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1  Note  to  the  Reader 


This  technical  report  is  the  extended  version  of  the  paper  “Polynomial-  Time  Verification  of 
PCTL  Properties  of  MDPs  with  Convex  Uncertainties”  presented  at  the  25th  International 
Conference  on  Computer  Aided  Verification,  CAV  2013.  The  report  extends  the  confer¬ 
ence  submission  with  the  following  additional  material,  which  was  not  added  due  to  space 
limitations: 

1.  Appendix  A:  Results  on  convex  optimization. 

2.  Appendix  B:  Alternative  verification  procedure  based  on  Value  Iteration. 

3.  Appendix  C:  Linear  programming  formulation  to  verify  the  property  including  the 
Until  operator  on  the  running  example  introduced  in  the  paper. 

4.  Appendix  D:  Case  study  on  the  verification  of  the  Dining  Philosopher  problem. 
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Abstract.  We  address  the  problem  of  verifying  Probabilistic  Computation  Tree 
Logic  (PCTL)  properties  of  Markov  Decision  Processes  (MDPs)  whose  state  transi¬ 
tion  probabilities  are  only  known  to  lie  within  uncertainty  sets.  We  first  introduce  the 
model  of  Convex-MDPs  (CMDPs),  i.e.,  MDPs  with  convex  uncertainty  sets.  CMDPs 
generalize  Interval-MDPs  (IMDPs)  by  allowing  also  more  expressive  (convex)  de¬ 
scriptions  of  uncertainty.  Using  results  on  strong  duality  for  convex  programs,  we 
then  present  a  PCTL  verification  algorithm  for  CMDPs,  and  prove  that  it  runs  in 
time  polynomial  in  the  size  of  a  CMDP  for  a  rich  subclass  of  convex  uncertainty 
models.  This  result  allows  us  to  lower  the  previously  known  algorithmic  complexity 
upper  bound  for  IMDPs  from  co-NP  to  PTIME.  Using  the  proposed  approach,  we 
verify  a  consensus  protocol  and  a  dynamic  configuration  protocol  for  IPv4  addresses. 


1  Introduction 

Stochastic  models  like  Discrete-Time  Markov  Chains  (DTMCs)  [1]  and  Markov  Decision 
Processes  (MDPs)  [2]  are  used  to  formally  represent  systems  that  exhibit  probabilistic  be¬ 
haviors.  These  systems  need  quantitative  analysis  [3]  to  answer  questions  such  as  “what  is  the 
probability  that  a  request  will  be  eventually  served?” .  Properties  of  these  systems  can  be  ex¬ 
pressed  and  analyzed  using  logics  such  as  Probabilistic  Computation  Tree  Logic  (PCTL)  [4] 
—  a  probabilistic  logic  derived  from  CTL  -  as  well  as  techniques  for  probabilistic  model 
checking  [5].  These  methods  often  rely  on  deriving  a  probabilistic  model  of  the  underly¬ 
ing  process,  hence  the  formal  guarantees  they  provide  are  only  as  good  as  the  estimated 
model.  In  a  real  setting,  these  estimations  are  affected  by  uncertainties  due,  for  example,  to 
measurement  errors  or  approximation  of  the  real  system  by  mathematical  models. 

Interval-valued  Discrete-Time  Markov  Chains  (IDTMCs)  have  been  introduced  to  cap¬ 
ture  modeling  uncertainties  [6].  IDTMCs  are  DTMC  models  where  each  transition  prob¬ 
ability  lies  within  a  compact  range.  Two  semantic  interpretations  have  been  proposed  for 
IDTMCs  [7]:  Uncertain  Markov  Chains  (UMCs)  and  Interval  Markov  Decision  Processes 
(IMDPs).  An  UMC  is  interpreted  as  a  family  of  DTMCs,  where  each  member  is  a  DTMC 
whose  transition  probabilities  lie  within  the  interval  range  given  in  the  UMC.  In  IMDPs, 
the  uncertainty  is  resolved  through  non-determinism.  Each  time  a  state  is  visited,  a  transi¬ 
tion  distribution  within  the  interval  is  adversarially  picked,  and  a  probabilistic  step  is  taken 
accordingly.  Thus,  IMDPs  model  a  non-deterministic  choice  made  from  a  set  of  (possibly 
uncountably  many)  choices.  In  this  paper  we  do  not  consider  UMCs  and  focus  on  IMDPs. 

An  upper-bound  on  the  complexity  of  model  checking  PCTL  properties  on  IMDPs  was 
previously  shown  to  be  co-NP  [8].  This  result  relies  on  the  construction  of  an  equivalent 
MDP  that  encodes  all  behaviors  of  the  IMDP.  For  each  state  in  the  new  MDP,  the  set  of 
transition  probabilities  is  equal  to  the  Basic  Feasible  Solutions  (BFS)  of  the  set  of  inequalities 
specifying  the  transition  probabilities  of  the  IMDP.  Since  the  number  of  BFS  is  exponential 


Table  1:  Known  Upper-Bound  on  the  Complexity  of  PCTL  Model  Checking. 


Model 

DTMC  [4] 

IMDP  [8] 

IMDP/CMDP  [ours] 

Complexity 

PTIME 

co-NP 

PTIME 

in  the  number  of  states  in  the  IMDP,  the  equivalent  MDP  can  have  size  exponential  in  the 
size  of  the  IMDP.  In  this  paper,  we  describe  a  polynomial-time  algorithm  (in  both  size  of  the 
model  and  size  of  the  formula)  based  on  Convex  Programming  (CP)  for  the  same  fragment 
of  PCTL  considered  in  [7,8]  (the  Bounded  Until  operator  is  disallowed).  This  shows  that 
the  problem  is  in  the  complexity  class  PTIME.  With  Bounded  Until ,  the  time  complexity  of 
our  algorithm  only  increases  to  pseudo-polynomial  in  the  maximum  integer  time  bound. 

An  interval  model  of  uncertainty  may  appear  to  be  the  most  intuitive.  However,  there 
are  significant  advantages  in  accommodating  also  more  expressive  (and  less  pessimistic) 
uncertainty  models.  In  [9],  a  financial  portfolio  optimization  case-study  is  analyzed  in  which 
uncertainty  arises  from  estimating  the  asset  return  rates.  The  authors  claim  that  the  interval 
model  is  too  conservative  in  this  scenario,  because  it  would  suggest  to  invest  the  whole 
capital  into  the  asset  with  the  smallest  worst-case  return.  The  ellipsoidal  model  proposed 
in  that  paper  returns  instead  the  more  profitable  strategy  of  spreading  the  capital  across 
multiple  assets.  Further,  depending  on  the  field,  researchers  use  different  models  to  represent 
uncertainty.  Maximum  likelihood  models  are  often  used,  for  example,  to  estimate  chemical 
reaction  parameters  [10].  To  increase  modeling  expressiveness,  we  introduce  the  model  of 
Convex-MDP  (CMDP),  i.e. ,  an  MDP  whose  state  transition  probabilities  are  only  known 
to  lie  within  convex  uncertainty  sets.  The  proposed  algorithms  can  be  extended  to  verify 
CMDPs  for  all  the  models  of  uncertainty  that  satisfy  a  technical  condition  introduced  later  in 
the  paper,  while  maintaining  the  same  complexity  results  proven  for  IMDPs.  This  condition 
is  not  a  limitation  in  practical  scenarios,  and  we  show  that  all  the  models  in  the  wide 
and  relevant  class  of  convex  uncertainty  sets  introduced  in  [11]  (e.g.  interval,  ellipsoidal 
and  likelihood  models)  satisfy  it.  Heterogeneous  models  of  uncertainty  can  then  be  used 
within  the  same  CMDP  to  represent  different  sources  of  uncertainty.  We  also  note  that  the 
complexity  results  presented  in  [7]  and  [8]  cannot  be  trivially  extended  to  verifying  CMDPs. 
This  is  because  BFS  are  not  defined  for  generic  convex  inequalities,  so  the  construction  of 
an  equivalent  MDP  would  not  be  possible.  The  complexity  results  are  compared  in  Table  1. 

To  summarize,  the  contributions  of  this  paper  are  as  follows. 

1.  We  give  a  polynomial-time  algorithm  for  model  checking  PCTL  properties  (without 
Bounded  Until)  on  IMDPs.  This  improves  the  co-NP  result  in  [8]  to  PTIME. 

2.  We  extend  the  algorithm  to  full  PCTL  (with  Bounded  Until)  and  show  that  its  time 
complexity  becomes  pseudo-polynomial  in  the  maximum  integer  bound  in  Bounded  Until. 

3.  We  show  that  our  complexity  results  extend  to  Convex-MDPs  (CMDPs)  for  a  wide  and 
expressive  subclass  of  the  convex  models  of  uncertainty. 

4.  We  demonstrate  the  relevance  of  our  approach  with  case  studies,  where  a  small  uncertainty 
in  the  probability  transitions  indeed  yields  a  significant  change  in  the  verification  results. 

The  paper  is  organized  as  follows.  Section  2  gives  background  on  MDPs,  PCTL,  and  the 
analyzed  uncertainty  models.  Section  3  presents  related  work  in  the  fields  of  verification  and 
control.  Section  4  gives  an  overview  of  the  proposed  approach.  In  Section  5,  we  describe  the 
proposed  algorithm  in  detail  and  prove  the  PTIME  complexity  result.  Section  6  describes 
two  case  studies,  and  we  conclude  and  discuss  future  directions  in  Section  7. 


2  Preliminaries 


Definition  2.1.  A  Probability  Distribution  (PD)  over  a  finite  set  Z  of  cardinality  n  is  a 
vector  n  £  ]R”  satisfying  0  <  n  <  1  and  1T  p,  =  1.  The  element  p[i]  represents  the  probability 
of  realization  of  the  event  Zi .  We  call  Dist(Z)  the  set  of  distributions  over  Z. 

2.1  Convex  Markov  Decision  Process  (CMDP) 

Definition  2.2.  A  CMDP  is  a  tuple  Me  =  {S,  So,  A,  17,  F ,  A,  X ,  L),  where  S  is  a  finite  set 
of  states  of  cardinality  N  =  |S|,  So  is  the  set  of  initial  states,  A  is  a  finite  set  of  actions 
(M  =  \A\),  12  is  a  finite  set  of  atomic  propositions,  F  is  a  finite  set  of  convex  sets  of 
transition  PDs,  A  :  S  — X  2A  is  a  function  that  maps  each  state  to  the  set  of  actions  available 
at  that  state,  1  =  S  x  1  4  J  is  a  function  that  associates  to  state  s  and  action  a  the 
corresponding  convex  set  Ff  £  F  of  transition  PDs,  and  L  :  S  — X  217  is  a  labeling  function. 

The  set  Ff  =  Dist).(S)  represents  the  uncertainty  in  defining  a  transition  distribution  for 
A Ac  given  state  s  and  action  a.  We  call  f“  £  Ff  an  observation  of  this  uncertainty.  Also, 
fa  g  jgjv  ancj  we  can  conect  vectors  f“,Vs  £  S  into  an  observed  transition  matrix  Fa  £ 
RNxJv.  Abusing  terminology,  we  call  Fa  the  uncertainty  set  of  the  transition  matrices,  and 
Fa  £  Fa.  Ff  is  interpreted  as  the  row  of  Fa  corresponding  to  state  s.  Finally,  ffs  =  f“  [j] 
is  the  observed  probability  of  transitioning  from  s*  to  Sj  when  action  a  is  selected. 

A  transition  between  state  s  to  state  s'  in  a  CMDP  occurs  in  three  steps.  First,  an  action 
a  £  A{s)  is  chosen.  The  selection  of  a  is  nondeterministic.  Secondly,  an  observed  PD  f“  £  Ff 
is  chosen.  The  selection  of  f“  models  uncertainty  in  the  transition.  Lastly,  a  successor  state 
s'  is  chosen  randomly,  according  to  the  transition  PD  f“. 

fs°s  fs\ 

A  path  7r  in  A 4c  is  a  finite  or  infinite  sequence  of  the  form  so  — - — >  s i - >,  •  •  • ,  where 

Si  £  S,  ai  £  M(si)  and  /“*s  >  0  Vi  >  0.  We  indicate  with  77/in  ( ninf )  the  set  of  all  finite 

(infinite)  paths  of  Me-  7r[i]  is  the  ith  state  along  the  path  and,  for  finite  paths,  last{ tv)  is 
the  last  state  visited  in  n  £  IIfin.  IIS  =  {n  \  7r[0]  =  s}  is  the  set  of  paths  starting  in  state  s. 
To  model  uncertainty  in  state  transitions,  we  make  the  following  assumptions: 

Assumption  2.1.  Fa  can  be  factored  as  the  Cartesian  product  of  its  rows,  i.e.,  its  rows  are 
uncorrelated.  Formally,  for  every  a  £  A,  Fa  =  Ff0  x  •  •  •  x  FfN_1 .  In  [11]  this  assumption  is 
referred  to  as  rectangular  uncertainty. 

Assumption  2.2.  If  the  pi'obability  of  a  transition  is  zero  (non-zero)  for  at  least  one  PD 
in  the  uncertainty  set,  then  it  is  zero  (non-zero)  for  all  PDs. 

Formally,  3f“  £  Ff  :  ffs,  =  (^)O  =>  Vf“  £  Ff  :  /“g,  =  (^)O.  The  assumption  guarantees 
the  correctness  of  the  preprocessing  verification  routines  used  later  in  the  paper,  which  rely 
on  reachability  of  the  states  of  the  MDP  underlying  graph. 

We  determine  the  size  1Z  of  the  CMDP  Ale  as  follows.  A 4c  has  N  states,  O(M)  actions 
per  state  and  0(N2)  transitions  for  each  action.  Let  D)  denote  the  number  of  constraints 
required  to  express  the  rectangular  uncertainty  set  Ff  (e.g.  D“  =  0{2N)  for  the  interval 
model,  to  express  the  upper  and  lower  bounds  of  the  transition  probabilities  from  state  s  to 
all  states  s'  £  S),  and  D  =  max  D[.  The  overall  size  of  Me  is  thus  1Z  =  0(N2M+NMD). 

s£  S \cl£  A 

In  order  to  analyze  quantitative  properties  of  CMDPs,  we  need  a  probability  space  over 
infinite  paths  [12].  However,  a  probability  space  can  only  be  constructed  once  nondetermin¬ 
ism  and  uncertainty  have  been  resolved.  We  call  each  possible  resolution  of  nondeterminism 
an  adversary,  which  chooses  an  action  in  each  state  of  Me- 


Definition  2.3.  Adversary.  A  randomized  adversary  for  A ic  is  a  function  a  =  II fin  xi-> 
[0,1],  with  J2a{1  ast(ir))  a(7Tia)  =  U  and  a  £  A(last(ir))  if  a(ir,a)  >  0.  We  call  Adv  the  set 
of  all  adversaries  a  of  AAc- 

Conversely,  we  call  a  nature  each  possible  resolution  of  uncertainty,  i.e.,  a  nature  chooses  a 
transition  PD  for  each  state  and  action  of  A ic- 


Definition  2.4.  Nature.  Given  action  a  £  A,  a  randomized  nature  is  the  function  r]a  : 
Ilfin  x  Dist.(S)  ->  [0,1]  with  >  ??a(7r,f“)  =  1,  and  fsa  G  F?ast^)  if  7/°(7r,  f“)  >  0.  We 

call  Nat  the  set  of  all  natures  rja  of  AAc- 


An  adversary  a  (nature  rja)  is  memoryless  if  it  depends  only  on  lasti^K).  Also,  a  (?ya)  is 
deterministic  if  a(7r,a)  =  1  for  some  a  G  A(last(n))  (r]a(ir,  /“)  =  1  for  some  f“  G  ^). 


2.2  Models  of  Uncertainty 

We  only  consider  CMDPs  whose  transition  PDs  lie  in  uncertainty  sets  that  satisfy  Assump¬ 
tion  5.1  (introduced  later  for  ease  of  presentation).  This  assumption  holds  for  all  the  uncer¬ 
tainty  models  analyzed  in  [11].  We  report  results  for  the  interval,  likelihood  and  ellipsoidal 
models.  A  more  thorough  derivation  is  available  in  Appendix  A. 


Interval  Model.  Intervals  commonly  describe  uncertainty  in  transition  matrices: 

Ta  =  {fa  g  RW  |  0  <  f“  <  f“  <  f“  <  1,  T'T“  =  1}  (1) 

where  f“,  f°  G  R"  are  the  element-wise  lower  and  upper  bounds  of  f.  This  model  is  suitable 
when  the  transition  matrix  components  are  individually  estimated  by  statistical  data. 

Likelihood  Model.  This  model  is  appropriate  when  the  transition  probabilities  are  deter¬ 
mined  experimentally.  The  transition  frequencies  associated  to  action  a  £  A  are  collected  in 
matrix  Ha.  Uncertainty  in  each  row  of  Ha  can  be  described  by  the  likelihood  region  [13]: 

Fas  =  {f“  G  I  f“  >  0,  lTfs°  =  1,  £s,  hass,  log (/“,,)  >  ft}  (2) 

where  /?“  <  /3“maa,  =  £s,  /i“s,  log(/i“s,)  represents  the  uncertainty  level.  Likelihood  regions 
are  less  conservative  uncertainty  representations  than  intervals,  which  arise  from  projections 
of  the  uncertainty  region  onto  each  row  component. 

Ellipsoidal  Model.  Ellipsoidal  models  can  be  seen  as  a  second-order  approximation  of  the 
likelihood  model  [11].  Formally: 

Ta  =  {fa  g  RN  j  fa  >  Q  ^fa  =  1;  \\Ra  (fa  _  ha}  ||2  <  1,  0}  (3) 

where  matrix  represents  an  ellipsoidal  approximation  of  the  likelihood  Region  (2). 

Remark  2.1.  Each  set  F'f  within  the  same  CMDP  can  be  expressed  with  a  different  un¬ 
certainty  model  to  represent  different  sources  of  uncertainty. 

To  illustrate  our  results,  we  will  use  the  IMDP  Afc  hr  Figure  1,  with  S  =  {so  •  •  •  S3}, 
So  =  {so},  A  =  {a,  6},  fl  =  {w,tf},  A  :  {s0,sr,s2}  -X  {a}  ;{s3}  ->  {a,  6},  L  :  {s0,s3}  — X 
it  ;{s2}  —X  oj.  The  uncertainty  intervals  are  shown  next  to  each  transition.  For  example, 
K ,  =  {f“o  G  I  [0, 0.6, 0.2, 0]  <  f“  <  [0, 0.8, 0.5, 0],  £s,6S  /“ ,  -  1}. 


Table  2:  PCTL  semantics  for  CMDP 


s 

h 

True 

s 

h 

w 

iff  to  G  L(s) 

s 

N 

iff  s  \f=  b 

s 

H 

<t>  1  A  02 

iff  s  j=  0i  A  s  |=  02 

s 

N 

P\ap  [b] 

iff  Prob  {{tv  £  IIs(a,r]a)  \  7 r  \=  b})  N  p 

Vcc  G  Adv  and  r]a  £  Nat 

7 r 

N 

Xcj) 

iff  7r[l]  ^  b 

7T 

N 

4 >1 

02  iff  37  <  7c  |  7r[z]  |=  cj>2  A  Vj  <  i  ir  [j]  \=  cj> l 

7T 

b 

4>1  Ucj)  2 

iff  3k  >  0  7T  [=  0i  U-k4>2 

2.3  Probabilistic  Computation  Tree  Logic  (PCTL) 

We  use  PCTL,  a  probabilistic  logic  derived  from  CTL  which  includes  a  probabilistic  operator 
P  [4],  to  express  properties  of  CMDPs.  The  syntax  of  this  logic  is  defined  as  follows: 

<j>  ::=  True  |  w  |  -\<j>  |  cj>i  A  02  |  P^p  [b]  state  formulas 

b  ::=  Xcj)  |  cj) i  U-k(j) 2  |  U(j) 2  path  formulas 

where  w  €  12  is  an  atomic  proposition,  Me  {<,  <,  >,  >},  p  £  [0, 1]  and  k  £  N. 

Path  formulas  b  use  the  Next  (X),  Bounded  Until  ( U-k )  and  Unbounded  Until  fJ)  oper¬ 
ators.  These  formulas  are  evaluated  over  paths  and  only  allowed  as  parameters  to  the  PMp  [0] 
operator.  The  size  Q  of  a  PCTL  formula  is  defined  as  the  number  of  Boolean  connectives 
plus  the  number  of  temporal  operators  in  the  formula.  For  the  Bounded  Until  operator,  we 
denote  separately  the  maximum  time  bound  that  appears  in  the  formula  as  kmax •  Proba¬ 
bilistic  statements  about  MDPs  typically  involve  universal  quantification  over  adversaries 
a  £  Adv.  With  uncertainties,  for  each  action  a  selected  by  adversary  a,  we  will  further  quan¬ 
tify  across  nature  rj a  £  Nat  to  compute  the  worst  case  condition  within  the  action  range  of 
r\a,  i.e.,  the  uncertainty  set  Ff.  We  define  Ps{a,  ??a)[b]  =  Prob  ({7r  G  PLs{a,  rja)  \  ir  }=  b})  the 
probability  of  taking  a  path  ir  £  77 s  that  satisfies  b  under  adversary  a  and  nature  7/°.  If  a 
and  rja  are  Markov  deterministic  in  state  s,  we  write  Ps(a,  f“),  where  a  and  f“  are  the  action 
and  resolution  of  uncertainty  that  are  deterministically  chosen  at  each  execution  step  by  a 
and  770.  P™ax[il>\  (p™» [0] )  denote  the  maximum  (minimum)  probability  Ps{a,  r/a)[ip]  across 
all  adversaries  a  £  Adv  and  natures  pa  £  Nat,  and  the  vectors  Pmaa:[7/)],  Pmm[-0]  £  R" 
collect  these  probabilities  Vs  £  S.  The  semantics  of  the  logic  is  reported  in  Table  2,  where 
we  write  \=  instead  of  \=  Adv, Nat  for  simplicity. 

For  ease  of  computation,  we  would  like  to  restrict  our  attention  to  memoryless  and  deter¬ 
ministic  adversaries  and  natures  to  compute  quantitative  probabilities,  i.e.,  solve  problems: 

Prxm=  max  max  Ps (a,  f“ ) [tp\  or  =  min  min  Ps(a,  f“)[V>]  (4) 

aeyt(s)f“eJr“  ae.4(s)f“e.F“ 

We  extend  a  result  from  [14]  to  prove  that  this  is  possible. 

Proposition  2.1.  Given  a  CMDP  M.c  and  a  target  state  St  G  S,  there  always  exist  de¬ 
terministic  and  memoryless  adversaries  and  natures  for  A ic  that  achieve  the  maximum 
(minimum)  probabilities  of  reaching  St,  if  A  is  finite  and  the  inner  optimization  in  Problem 
(4)  always  attains  its  optimum  cr*(a)  over  the  sets  P“,Vs  G  S,\/a  G  M(s),  i.e.,  there  exists 
a  finite  feasible  if  £  Ff  such  that  Ps(a,  f“)[V>]  =  <7* (a). 


Sketch  of  proof .  The  proof  is  divided  into  two  parts.  First,  we  prove  the  existence  of  an 
adversary  and  a  nature  that  achieve  the  maximum  (minimum)  probabilities  of  reaching  ts, 
using  Banach  fixed-point  theorem  [14].  Second,  we  prove  that  at  least  one  such  adversary 
and  nature  is  memoryless  and  deterministic.  The  proof  extends  the  one  in  Puterman  [14], 
Theorem  6.2.10.  We  need  to  prove  that  Problem  (4)  always  attains  the  maximum  (minimum) 
over  the  feasibility  sets  i.e.,  Vs  €  S,\/a  €  M(s),3f“  €  Tf  :  |  |f“  1 12  <  oo,  Ps(a,  ff)[ip]  = 
o'* (a).  This  is  indeed  true  for  all  the  convex  sets  Ff  considered  in  this  paper.  The  interval 
and  ellipsoidal  models  result  in  compact  sets  F'f  on  which  Weierstrass  theorem  holds.  For 
the  likelihood  model  we  use  the  notion  of  consistency,  which  guarantees  the  existence  and 
uniqueness  of  a  point  in  T'f  where  the  optimum  is  attained.  □ 

The  verification  of  a  PCTL  state  formula  <p  can  be  viewed  as  a  decision  problem.  The 
verification  algorithm  V  needs  to  determine  whether  a  state  s  £  Sq  is  (or  is  not)  contained 
in  the  set  Sat(cp)  =  {s  G  S'  |  s  |=  cp}.  We  can  thus  define  the  following  properties  for  V: 

Definition  2.5.  Soundness  (^Completeness).  Algorithm  V  is  sound  (complete)  if: 

s  £  Satv(4>)  =>  s  £  Sat(cp)  (s  fL  Satv{4>)  =>  s  $  Sat(<p)) 
where  Satv(4>)  is  the  computed  satisfaction  set,  while  Sat(<p)  is  the  actual  satisfaction  set. 

Algorithms  to  verify  non-probabilistic  formulas  are  sound  and  complete,  because  they  are 
based  on  reachability  analysis  over  the  finite  number  of  states  of  A ic  [15].  Conversely,  we 
will  show  in  Section  5  that  algorithms  to  verify  probabilistic  formulas  <p  =  [ip]  in  the 
presence  of  uncertainties  require  to  solve  convex  optimization  problems  over  the  set  ffi.  of  the 
real  numbers.  Optima  of  these  problems  can  be  arbitrary  real  numbers,  so,  in  general,  they 
can  be  computed  only  to  within  a  desired  accuracy  e^.  We  consider  an  algorithm  to  be  sound 
and  complete  if  the  error  in  determining  the  satisfaction  probabilities  of  <p  is  bounded  by 
such  a  parameter  e<j,  since  the  returned  result  will  still  be  accurate  enough  in  most  settings. 


3  Related  Work 

Probabilistic  model  checking  tools  such  as  PRISM  [5]  have  been  used  to  analyze  a  multitude 
of  applications,  from  communication  protocols  and  biological  pathways  to  security  problems. 
In  this  paper,  we  further  consider  uncertainties  in  the  probabilistic  transitions  of  the  MDP 
for  model  checking  PCTL  specifications.  Prior  work  [6-8, 16]  in  similar  verification  problems 
also  dealt  with  uncertainties  in  the  probabilistic  transitions.  However,  they  considered  only 
interval  models  of  uncertainty,  while  we  incorporate  more  expressive  models  such  as  ellip¬ 
soidal  and  likelihood.  Further,  we  consider  nature  as  adversarial  and  study  how  it  affects  the 
MDP  execution  in  the  worst  case.  The  developers  of  PARAM  [17]  consider  instead  uncer¬ 
tainties  as  possible  values  that  parameters  in  the  model  can  take,  and  synthesize  the  optimal 
parameter  values  to  maximize  the  satisfaction  probability  of  a  given  PCTL  specification. 

We  improve  the  previously  best-known  complexity  result  of  co-NP  in  [8]  to  PTIME, 
for  the  fragment  of  PCTL  without  U-k .  For  the  full  PCTL  syntax,  our  algorithm  runs  in 
0(poly( 71)  x  Q  x  kmax)  time,  where  kmax  is  the  maximum  bound  in  U-k.  This  result  is 
pseudo-polynomial  in  kmax,  i.e.,  polynomial  (exponential)  if  kmax  is  counted  in  its  unary 
(binary)  representation.  Conversely,  classical  PCTL  model  checking  for  DTMCs  [4]  runs  in 
time  polynomial  in  kmax  counted  in  its  binary  representation.  The  difference  stems  from 
the  computation  of  the  set  Sat  (PNp  [(pi  U-k<p2\ ) .  For  (certain)  MDPs,  this  computation 
involves  raising  the  transition  matrices  Fa,\/a  £  A  to  the  kth  power,  to  model  the  evolution 
of  the  system  in  k  steps.  With  uncertainties,  we  cannot  do  matrix  exponentiation,  because 


F°  €  Ta  might  change  at  each  step.  However,  both  Q  and  kmax  are  typically  small  in 
practical  applications  [18, 19],  so  the  dominant  factor  for  runtime  is  the  size  of  the  model  7Z. 
We  note  that  the  complexity  results  of  [7]  and  [8]  can  be  extended  to  the  PCTL  with  U-k. 

The  convex  uncertainty  models  [11]  analyzed  in  this  paper  have  been  considered  recently 
in  the  robust  control  literature.  In  [20],  an  algorithm  is  given  to  synthesize  a  robust  optimal 
controller  for  an  MDP  to  satisfy  a  Linear  Temporal  Logic  (LTL)  specification  where  only 
one  probabilistic  operator  is  allowed.  Their  technique  first  converts  the  LTL  specification  to 
a  Rabin  automaton  (which  is  worst-case  doubly  exponential  in  the  size  of  the  LTL  formula), 
and  composes  it  with  the  MDP.  Robust  dynamic  programming  is  then  used  to  solve  for  the 
optimal  control  policy.  We  consider  PCTL,  which  allows  nested  probability  operators,  and 
propose  an  algorithm  which  is  polynomial  both  in  the  size  of  the  model  and  of  the  formula. 

In  [21],  the  robustness  of  PCTL  model  checking  is  analyzed  based  on  the  notion  of  an 
Approximate  Probabilistic  Bisimulation  (APB)  tailored  to  the  finite-precision  approximation 
of  a  numerical  model.  We  instead  verify  MDPs  whose  transition  probabilities  are  affected 
by  uncertainties  due  to  estimation  errors  or  imperfect  information  about  the  environment. 


4  Probabilistic  Model  Checking  with  Uncertainties 

We  define  the  problem  under  analysis,  and  overview  the  proposed  approach  to  solve  it. 

PCTL  model  checking  with  uncertainties.  Given  a  Markov  Decision  Process  model  with 
convex  uncertainties  Ate  of  size  7Z  and  a  PCTL  formula  (j)  of  size  Q  over  a  set  of  atomic 
propositions  C,  verify  cj>  over  the  uncertainty  sets  J“  £  J  of  Ate- 

As  in  verification  of  CTL  [22],  the  algorithm  traverses  bottom- up  the  parse  tree  for  (f), 
recursively  computing  the  set  Sat(cj)')  of  states  satisfying  each  sub- formula  (j)' .  At  the  end 
of  the  traversal,  the  algorithm  computes  the  set  of  states  satisfying  (f>  and  it  determines  if 
s  \=  (j)  by  checking  if  s  G  Sat  (4>).  For  the  non-probabilistic  PCTL  operators,  the  satisfying 
states  are  computed  as:  Sat  (True)  =  S ,  Sat(oj )  =  {s  €  S  \  u  €  L(s)}1  Sat(-«j))  =  S\Sat((f) 
and  Sat((j>i  A  fo)  =  Sat(cj>  1)  (~l  Sat((f >2).  For  the  probabilistic  operator  P  N  [if],  we  compute: 

Sat  ( P,p  M)  =  {s  G  S  |  PT* (</>)  <p}  Sat  ( P>p  [if])  =  {s  G  S'  |  P™n  W  >p}  (5) 

In  this  paper,  we  propose  polynomial-time  routines  to  compute  Sets  5  for  MDPs  whose 
transition  matrices  Fa  are  only  known  to  lie  within  convex  uncertainty  sets  Fa,  Va  €  A. 

Using  Proposition  2.1,  the  proposed  routines  encode  the  transitions  of  Ale  under  the 
sets  of  deterministic  and  nrenroryless  adversaries  and  natures  into  convex  programs  and 
solve  them.  From  the  returned  solution,  it  is  then  possible  to  determine  the  quantitative 
satisfaction  probabilities  Pfnax[if]  (or  Pfam[if])  Vs  £  S,  which  get  compared  in  linear  time 
to  the  threshold  p  to  compute  the  set  Sat  (PMp  [ip])-  To  prove  the  polynomial-time  complexity 
of  the  model-checking  algorithm,  we  use  the  following  key  result  from  convex  theory  [23] . 

Proposition  4.1.  Given  the  convex  program: 
min  /0(x) 

X 

s.t.  fi(x)  <0  *  =  1,  •  •  •  , m 

with  x  £  Kn  and  fi,i  =  0,  •  •  •  ,m  convex  functions,  the  optimum  <j*  can  he  found  to  within 
±ed  in  time  complexity  that  is  polynomial  in  the  size  of  the  problem  ( n,m )  and  log(l/ed). 


We  are  now  ready  to  state  the  main  contribution  of  this  paper: 


Theorem  4.1.  Complexity  of  PCTL  model-checking  for  CMDPs. 

1.  The  problem  of  verifying  if  a  CMDP  A ic  of  size  1Z  satisfies  a  PCTL  formula  (j>  without 
U^k  is  in  PTIME. 

2.  A  formula  </>'  with  U-k  can  be  verified  with  time  complexity  O  (poly (TZ)  x  Qf  x  kmax), 
i.e.,  pseudo-polynomial  in  the  maximum  time  bound  kmax  ofU-k. 

Sketch  of  proof.  The  proof  is  constructive.  Our  verification  algorithm  parses  (f>  in  time  linear 
in  the  size  Q  of  </>  [22],  computing  the  satisfiability  set  of  each  operator  in  4>.  For  the  non- 
probabilistic  operators,  satisfiability  sets  can  be  computed  in  time  polynomial  in  7 Z  using 
set  operations,  i.e.,  set  inclusion,  complementation  and  intersection.  For  the  probabilistic 
operator,  we  leverage  Proposition  4.1  and  prove  that  the  proposed  verification  routines:  1) 
solve  a  number  of  convex  problems  polynomial  in  TZ;  2)  generate  these  convex  programs  in 
time  polynomial  in  TZ.  The  correctness  and  time-complexity  for  formulas  involving  Next  and 
Unbounded  Until  operators  are  formalized  in  Lemma  5.1  and  5.2  (Section  5.1  and  5.2).  It 
thus  follows  that  the  overall  algorithm  runs  in  time  polynomial  in  TZ  and  in  the  size  of  f>. 
Finally,  Lemma  5.3  formalizes  the  results  related  to  the  Bounded  Until  operator.  □ 

5  Verification  Routines 

We  detail  the  routines  used  to  verify  the  probabilistic  operator  P.  We  only  consider  prop¬ 
erties  in  the  form  =  P<p[ip\,  but  the  results  can  trivially  be  extended  to  (j)  =  P^p[^]  by 
replacing  “max”  with  “min”  in  the  optimization  problems  below. 

5.1  Next  Operator 

We  verify  property  =  Pmp[X4>i]  on  a  CMDP  of  size  TZ.  First,  the  set  Syes  =  Sat((j> i)  is 
computed.  Second,  for  all  state  s  G  S,  the  algorithm  evaluates  Equation  (4)  by  solving: 

P™ax[X<j) i]=  max  maxEs, eS»„ (6) 

a6-4(s)rf 

The  inner  max  is  a  convex  program  since  Ff  is  convex.  The  sets  Ff  can  be  expressed 
with  heterogeneous  uncertainty  models,  since  each  problem  is  independent  from  the  others. 
Finally,  the  computed  probabilities  are  compared  to  p  to  select  the  states  that  satisfy  (f>. 

Lemma  5.1.  The  routine  to  verify  the  Next  operator  is  sound,  complete  and  guaranteed  to 
terminate  with  algorithmic  complexity  that  is  polynomial  in  the  size  TZ  of  Ale- 

Proof.  Problem  (6)  has  one  “inner”  convex  program  Vs  €  S  and  Va  £  M(s),  for  a  total 
of  0(MN)  problems.  Each  problem  has  O(N)  unknowns,  representing  the  probability  of 
transitioning  from  state  s  to  state  s'  for  s'  €  Syes.  It  has  0(N  + 1)  constraints  to  guarantee 
that  the  solution  lies  in  the  probability  simplex,  and  D“  constraints  to  enforce  the  solution 
to  be  within  the  uncertainty  set  Ff .  The  total  number  of  constraints  is  thus  linear  in  TZ. 
Using  Proposition  4.1,  each  inner  problem  is  solved  in  time  polynomial  in  TZ.  Soundness  and 
completeness  also  follow  directly  from  Proposition  4.1,  which  states  that  the  optimum  of 
Problem  (6)  can  be  found  to  within  the  desired  accuracy  in  time  linear  in  log(l/ed).  □ 

We  verify  <f>  =  P<o.4[Tu;]  in  the  example  in  Figure  1.  Trivially,  Syes  =  {s2j--  We  solve 
Problem  (6)  for  all  a  €  A  and  s  €  S.  As  an  example,  for  state  s0  and  action  a,  we  solve: 

j-ta.max  n 

Pso  =  max  f 02 

701 >/02 

s.t.  0.6  <  foi  <  0.8;  0.2  <  f02  <  0.5;  f0 1  +  /02  =  1 
and  get  Pf’max[Xco\  =  0.4.  Overall,  we  get  Pmaa:[Tw]=[0.4,  0.5, 0, 0.6],  so  Sat(<j))={s0,  s2}- 


5.2  Unbounded  Until  Operator 


We  verify  <f>  =  P<p[<j)ihl<f>2\  on  a  CMDP  of  size  7Z.  First,  the  sets  Syes  =  Sat  (P>i  [<j>iU  fa}) , 
Sno  =  Sat  (P<o[fal(fa])  and  S1  =  S  \  (. Sno  U  Syes)  are  precomputed  in  time  polynomial  in 
TZ  using  conventional  reachability  routines  over  the  CMDP  underlying  graph  [15].  Second, 
Equation  (4)  is  evaluated  for  all  s  €  S  at  the  same  time  using  the  Convex  Programming  pro¬ 
cedure  described  next.  Finally,  the  computed  probabilities  are  compared  top.  An  alternative 
verification  procedure,  based  on  Value  Iteration,  can  be  found  in  Appendix  B. 


Convex  Programming  Procedure  (CP).  We  start  from  the  classical  LP  formulation  to 
solve  the  problem  without  the  presence  of  uncertainty  [15]: 

min  xTl 

X 

s.t.  xs  =0;  xs  =  1;  Vs  €  Sno;  s  £  Syes;  (7) 

xs  >  xTff  Vs  £  S7 ,Vo  £  A(s) 

where  Pmax[faL(fa\  =  x*  is  computed  solving  only  one  LP.  Problem  (7)  has  N  unknowns 
and  N  —  Q  +  MQ  constraints,  where  Q  =  IS1’  |  =  O(N),  so  its  size  is  polynomial  in  TZ. 
Proposition  2.1  allows  us  to  rewrite  Problem  (7)  in  the  uncertain  scenario  as: 

min  xT  1 

X 

s.t.  xs  =0;  xa  =  1;  Vs  £  Sno;  Vs  £  Syes;  (8) 

xs  >  max  (xTfsa)  Vs  £  S'7,  Va  £  A(s) 

i.e. ,  we  maximize  the  lower  bound  on  xs  across  the  nature  action  range.  The  decision  variable 
of  the  inner  problem  is  f“  and  its  optimal  value  <r*(x)  is  parameterized  in  the  outer  problem 
decision  variable  x.  Problem  (8)  can  be  written  in  convex  form  for  an  arbitrary  uncertainty 
model  by  replacing  the  last  constraint  with  one  constraint  for  each  point  in  J~" .  However, 
this  approach  results  in  infinite  constraints  if  the  set  J-"“  contains  infinitely  many  points,  as 
in  the  cases  considered  in  the  paper.  We  solve  this  difficulty  using  duality,  which  allows  to 
rewrite  Problem  (8)  with  a  number  of  constraints  polynomial  in  7 Z.  We  start  by  replacing 
the  primal  inner  problem  in  the  outer  Problem  (8)  with  its  dual  Vs  £  S?  and  Va  £  _4(s): 

cr“(x)  =  ^max  xTf“  =>■  d“(x)  =  min  g( Ag,x)  (9) 

where  A®  is  the  (vector)  Lagrange  multiplier  and  is  the  feasibility  set  of  the  dual.  In 

the  dual,  the  decision  variable  is  A®  and  its  optimal  value  d“(x)  is  parameterized  in  x. 

The  dual  function  g{ A“,x)  and  the  set  are  convex  by  construction  in  A“  for  arbitrary 
uncertainty  models,  so  the  dual  is  convex.  Further,  since  also  the  primal  is  convex,  strong 
duality  holds,  i.e.,  cr“  =  d“,  Vx  £  because  the  primal  satisfies  Slater’s  condition  [24] 
for  any  non-trivial  uncertainty  set  ,F“.  Any  dual  solution  overestimates  the  primal  solution. 
When  substituting  the  primals  with  the  duals  in  Problem  (8),  we  drop  the  inner  optimization 
operators  because  the  outer  optimization  operator  will  find  the  least  overestimates,  i.e.,  the 
dual  solutions  Vs  £  S,  a  £  *4(s),  to  minimize  its  cost  function.  We  get  the  CP  formulation: 


min  xT  1 

min  xTl 

X 

x,  A 

s.t.  xs  =0;  Xs  =  1; 

s.t.  Xs  =0;  xs  =  1; 

Vs  £  Sno;Vs  £  Syes-, 

(10a) 

Xs  >  min  g  (A“,  x) 

Xs  >  g{ A“,x); 

Vs  £  S'7,  Va  £  Al(s); 

(10b) 

A:  £  Vas 

Vs  £  S7,  Va  £  Al(s) 

(10c) 

The  decision  variables  of  Problem  (10)  are  both  x  and  A®,  so  the  CP  formulation  is  convex 
only  if  the  dual  function  g( A®,  x)  is  jointly  convex  in  A®  and  x.  While  this  condition  cannot 
be  guaranteed  for  arbitrary  uncertainty  models,  we  prove  constructively  that  it  holds  for  the 
ones  considered  in  the  paper.  For  example,  for  the  interval  model,  Problem  (10)  reads: 


min  xTl 

x.Af 

s.t.  xs  =0;  xs  =  1; 

Vs  €  S"°;Vs  G  Syes 

>  A?,.  -  (fDTAS,s  +  (fI)TAt,s; 

Vs  €  S? ,  Va  G  A(s); 

x  +  A2,s  —  A3  s  —  Ai  sl  =  0; 

Vs  G  S'7,  Va  G  Vl(s); 

A5,s  >  0,  A|,a  >  0 

Vs  G  S?,Va  G  A(s) 

which  is  an  LP,  so  trivially  jointly  convex  in  x  and  A®.  Analogously,  Problem  (10)  for  the 
ellipsoidal  model  is  a  Second-Order  Cone  Program  (SOCP),  as  reported  in  Appendix  A,  so 
again  jointly  convex  in  x  and  A®.  For  the  likelihood  model,  Problem  (10)  reads: 


min  xJT 

zs,A“ 


s.t.  xs  =0;  xs  =  1; 

Vs  G  S"°;Vs  G  Syes; 

(Ha) 

Xs  >  Ax^  (1  +  j3a  )A2iS+A2,s  hss’  l°g  (  AJ  s~xs,  )  ’ 

Vs  G  S?,  Va  G  Vl(s); 

(lib) 

Ai,s  >  max  xs' ;  A2,s  >  0 

Vs  G  S?,  Va  G  Vl(s) 

(11c) 

We  prove  its  joint  convexity  in  x  and  A®  as  follows.  The  cost  function  and  Constraints  (11a)- 
(11c)  are  trivially  convex.  Constraint  (lib)  is  generated  by  a  primal-dual  transformation, 
so,  according  to  convex  theory,  it  is  convex  in  the  dual  variables  A®  by  construction.  Convex 
theory  also  guarantees  that  the  affine  subtraction  of  x  from  A®  s  preserves  convexity,  given 
A“s  >  max  av,Vs  G  S  in  Constraint  (11c),  so  we  conclude  that  Problem  (11)  is  convex. 

For  general  CMDPs,  we  will  assume: 

Assumption  5.1.  Given  a  CMDP  A 4c,  for  all  convex  uncertainty  sets  Ff  G  F,  the  dual 
function  g( A“,x)  in  Problem  (9)  is  jointly  convex  in  both  A®  and  x. 

According  to  Proposition  4.1,  Problem  (10)  can  thus  be  solved  in  polynomial  time. 
Also  for  this  formulation,  ~Pmax\(j)iU(f>2 ]  =  x*,  so  all  the  satisfaction  probabilities  can  be 
computed  by  solving  only  one  convex  problem.  Finally,  we  note  that  we  can  combine  models 
of  uncertainty  different  from  one  another  within  a  single  CP  formulation,  since  each  dual 
problem  is  independent  from  the  others  according  to  Assumption  2.1.  As  an  example,  if  both 
the  interval  and  ellipsoidal  models  are  used,  the  overall  CP  formulation  is  an  SOCP. 

Lemma  5.2.  The  routine  to  verify  the  Unbounded  Until  operator  is  sound,  complete  and 
guaranteed  to  terminate  with  algorithmic  complexity  polynomial  in  the  size  7Z  of  Me,  if  Me 
satisfies  Assumption  5.1. 

Proof.  The  routine  solves  only  one  convex  program,  generated  in  time  polynomial  in  TZ  as 
follows.  We  formulate  Constraints  (10b)  and  (10c)  Vs  €  S?  and  a  €  -4(s),  i.e.,  O(MQ) 
constraints,  where  Q  =  |S',|  =  O(N).  They  are  derived  from  MQ  primal-dual  transfor¬ 
mations  as  in  Equation  (9).  Each  primal  problem  has  N  unknowns,  N  +  1  constraints  to 
represent  the  probability  simplex  and  Z?“  constraints  to  represent  the  uncertainty  set  Ff. 
From  duality  theory,  the  corresponding  dual  inner  problem  has  N  +  1  +  D®  unknowns  and 
2 N  +  1  +  -D®  constraints.  Overall,  Problem  (10)  has  O  (( N  +  1  +  D)MQ)  more  unknowns 
and  O  ((2 N  +  1  +  D)MQ)  more  constraints  of  Problem  (7),  so  its  size  is  polynomial  in  TZ. 


If  A ic  satisfies  Assumption  5.1,  Problem  (10)  is  convex.  Using  Proposition  4.1,  we  conclude 
that  it  can  be  solved  in  time  polynomial  in  1Z.  Finally,  when  strong  duality  holds  for  the 
transformation  in  Equation  (9),  soundness  and  completeness  of  the  final  solution  are  pre¬ 
served  because  the  dual  and  primal  optimal  value  of  each  inner  problem  are  equivalent.  □ 

We  verify  <f>  =  P>o.3[  d  U  u>  ]  on  the  example  in  Figure  1.  Problem  (10)  written  with  the 
data  of  the  model  has  19  variables  and  11  constraints,  and  it  can  be  found  in  Appendix  C. 
The  solution  reads:  Pmra[  iWw]  =  [0.2,  0, 1, 0.32],  and,  in  conclusion,  Sat((j> )  =  {s2>  s3}. 


5.3  Bounded  Until  Operator 

We  present  the  routine  to  verify  property  <f  =  PMp[</> iU-k(j>2\  on  a  CMDP  of  size  1Z.  First, 

the  set  Syes  d=  Sat{<t>2 ),  Sno  d=  S  \  (Sat(<j> i)  U  Sat(fa))  and  S?  =  S  \  (Sno  U  Syes)  are 
precomputed.  Second,  the  maximum  probabilities  Pmax  [^]  =  xfc  =  Gk( 0)  to  satisfy  <f>  are 
computed  for  all  states  s  £  S  applying  k  times  mapping  G  defined  as: 


i  ii  /  i  —  1  \ 

x  =  G  (x  ) 


0;  1;  Vs  €  Sno;  Vs  €  Syes; 
0;  Vs  £  S7  A  i  =  0; 
max  max  (xI_1)Tf“  Vs  £  S?  A  i  >  0 

a6^(s)f“6JrJ 


(12) 


and  x  1  =  0  €  RN.  Finally,  the  computed  probabilities  are  compared  to  the  threshold  p. 


Lemma  5.3.  The  routine  to  verify  the  Bounded  Until  operator  is  sound,  complete  and 
guaranteed  to  terminate  with  algorithmic  complexity  that  is  polynomial  in  the  size  1Z  of  Me 
and  pseudo-polynomial  in  the  time  bound  k  ofU-k. 


Proof.  The  proof  of  polynomial  complexity  in  7 Z  is  similar  to  the  one  for  the  Next  Operator. 
Further,  the  pseudo-polynomial  complexity  in  k  comes  from  applying  Mapping  (12)  k  times. 
While  each  inner  problem  is  solved  with  accuracy  ±£jran  in  time  linear  in  log(l/ejnn)  by 
Proposition  4.1,  we  also  need  to  prove  the  soundness  and  completeness  of  the  overall  solution, 
since  the  firm-approximations  in  x4,  Vi,  get  propagated  at  each  iteration  and  the  error  might 
get  amplified.  We  call  e4  the  error  accumulated  at  step  i  for  state  s,  x\  =  x4  id  +  e4 ,  where 
x4  id  is  the  solution  with  no  error  propagation,  and  ek  the  error  in  the  final  solution.  Also, 

f“’4  €  J7"  is  the  optimal  solution  of  the  inner  problem  at  step  i.  We  solve  this  difficulty  by 
noting  that  the  optimal  value  of  the  inner  problem  is  computed  by  multiplying  vector  x4  by 
f“’4  £  Pf,  with  lTf“  =  l,Vf“  £  Pf,\/a  £  A(s).  At  the  first,  second  and  ith  iteration: 


1  _  1  .  1  _  ca.  1  0  . 

xs  ^s,id  i  €s  Is  X  Cjnn 

2  -  na, 2  1  .  -  nd,2  (nCl,l  0  .  1  'j  f 

*^s  —  *-s  i  ^xnn  —  ^-s  \  As  I  ^inn  )  i  c 


fa, 2nd, l  0  |  q 

s  Is  X  ~r  -"Cinn 

fa.i/na.i  —  1  _z  — 1  ,  /  •  -i  \  i\  i  na.ina.i—  1  na.l  _0  ,  • 

s  (is  ^  “1“  € inn  Is  Is  •  •  •  Is  H-  ^ i : 


T*  —  f°,i  „i-l  i  ,, 

s  As  i  cmn 


so  e4  increases  linearly  with  i.  The  desired  accuracy  td  of  the  final  solution  can  thus  be 
guaranteed  by  solving  each  inner  problem  with  accuracy  =  ed/k.  □ 

We  verify  (f)  =  P<o.6[d  U -1w]  in  the  example  in  Fig.  1.  Syes  =  {S2},  Sno  =  {si}.  Applying 
once  Mapping  (12),  we  get  Pmax[dU-1ijj]  =  [0.4, 0, 1,  0.6]  and  Sat((j>)  =  {so,  si,  53}- 


6  Case  Studies 

We  implemented  the  proposed  verification  algorithm  in  Python,  and  interfaced  it  with 
PRISM  [5]  to  extract  information  about  the  CMDP  model.  We  used  MOSEK  [25]  to  solve 


the  LPs  generated  for  the  interval  model  and  implemented  customized  numerical  solvers  for 
the  other  models  of  uncertainty.  The  implemented  tool  is  available  at  [26].  The  algorithm 
was  tested  on  all  the  case  studies  collected  in  the  PRISM  benchmark  suite  [27] .  Due  to  space 
limits,  we  report  two  of  them:  the  verification  of  a  consensus  protocol  and  of  a  dynamic  con¬ 
figuration  protocol  for  IPv4  addresses.  Further,  the  Dining  Philosopher  problem  is  verified 
in  Appendix  D.  The  goals  of  these  experiments  are  two-fold:  1)  quantitatively  evaluate  the 
impact  of  uncertainty  on  the  results  of  verification  of  PCTL  properties  of  CMDPs;  2)  assess 
the  scalability  of  the  proposed  approach  to  increasing  problem  size.  The  runtime  data  were 
obtained  on  a  2.4  GHz  Intel  Xeon  with  32GB  of  RAM. 

6.1  Consensus  Protocol 

Consensus  problems  arise  in  many  distributed  environments,  where  a  group  of  distributed 
processes  attempt  to  reach  an  agreement  about  a  decision  to  take  by  accessing  some  shared 
entity.  A  consensus  protocol  ensures  that  the  processes  will  eventually  terminate  and  take 
the  same  decision,  even  if  they  start  with  initial  guesses  that  might  differ  from  one  another. 

We  analyze  the  randomized  consensus  protocol  presented  in  [18,28].  The  protocol  guar¬ 
antees  that  the  processes  return  a  preference  value  v  £  {1,2},  with  probability  parameterized 
by  a  process  independent  value  R  (R  >  2)  and  the  number  of  processes  P.  The  processes 
communicate  with  one  another  by  accessing  a  shared  counter  of  value  c.  The  protocol  pro¬ 
ceeds  in  rounds.  At  each  round,  a  process  flips  a  local  coin,  increments  or  decrements  the 
shared  counter  depending  on  the  outcome  and  then  reads  its  value  c.  If  c  >  PR  (c  <  —PR), 
it  chooses  v  =  1  (v  =  2).  Note  that  the  larger  the  value  of  R ,  the  longer  it  takes  on  average 
for  the  processes  to  reach  the  decision.  Nondeterminism  is  used  to  model  the  asynchronous 
access  of  the  processes  to  the  shared  counter,  so  the  overall  protocol  is  modeled  as  an  MDP. 

We  verify  the  property  Agreement:  all  processes  must  agree  on  the  same  decision, 
i.e. ,  choose  a  value  v  £  { 1, 2}.  We  compute  the  minimum  probability  of  Agreement  and 
compare  it  against  the  theoretical  lower  bound  ( R  —  l)/2 R  [18].  In  PCTL  syntax: 

P™m  [4>\  :=  P™*n  (F  ({finished}  A  {all -Coins -equal A}))  (13) 

We  consider  the  case  where  one  of  the  processes  is  unreliable  or  adversarial,  i.e.,  it  throws 
a  biased  coin  instead  of  a  fair  coin.  Specifically,  the  probability  of  either  outcome  lies  in  the 
uncertainty  interval  [(1  —  u)po,  (1  +  u)po\,  where  po  =  0.5  according  to  the  protocol.  This 
setting  is  relevant  to  analyze  the  protocol  robustness  when  a  process  acts  erroneously  due 
to  a  failure  or  a  security  breach.  In  particular,  our  approach  allows  to  study  attacks  that 
deliberately  hide  under  the  noise  threshold  of  the  protocol.  In  such  attacks,  the  compromised 
node  defers  agreement  by  producing  outputs  whose  statistical  properties  are  within  the  noise 
tolerance  of  an  uncompromised  node,  so  that  it  is  harder  to  detect  its  malicious  behavior. 

Figure  2  shows  the  effect  of  different  levels  of  uncertainty  on  the  computed  probabilities 
for  P  =  4.  With  no  uncertainty  (u  =  0),  P™”  increases  as  R  increases,  because  a  larger  R 
drives  the  decision  regions  further  apart,  making  it  more  difficult  for  the  processes  to  decide 
on  different  values  of  v .  As  R  goes  to  infinity,  P™tn  approaches  the  theoretical  lower  bound 
limn^00(R  —  l)/2 R  =  0.5.  However,  even  with  a  small  uncertainty  (u  =  0.01),  P"im  soon 
decreases  for  increasing  R.  With  a  large  uncertainty  ( u  =  0.15),  P™m  quickly  goes  to  0.  A 
possible  explanation  is  that  the  faulty  process  has  more  opportunities  to  deter  agreement  for 
a  high  R,  since  R  also  determines  the  expected  time  to  termination.  Results  thus  show  that 
the  protocol  is  vulnerable  to  uncertainties.  This  fact  may  have  serious  security  implication, 
i.e.,  a  denial-of-service  attack  could  reduce  the  availability  of  the  distributed  service,  since 
a  compromised  process  may  substantially  alter  the  expected  probability  of  agreement. 


Fig.  2:  Value  of  Eq.  13  in  function  of  R  Fig.  3:  Scalability  of  the  CP  procedure, 

while  varying  the  uncertainty  level  u. 

Lastly,  we  study  the  scalability  of  the  CP  procedure,  by  evaluating  Equation  (13)  while 
sweeping  R  both  for  P  =  2  and  P  =  4.  We  use  MOSEK  [25]  to  solve  Problem  (10)  and  set 
the  Time  Out  (TO)  to  one  hour.  In  Figure  3,  we  plot  the  sum  (IV  +  T)  of  the  number  of 
states  (V)  and  transitions  (T)  of  the  CMDP,  which  are  independent  of  the  uncertainty  in 
the  transition  probabilities,  to  represent  the  model  size  (top),  the  sum  (V +C)  of  the  number 
of  variables  (V)  and  constraints  (C)  of  the  generated  LP  instances  of  Problem  (10)  (center), 
and  the  running  time  tc p  (bottom) .  V  +  C  always  scales  linearly  with  N  +  T  (the  lines  have 
the  same  slope),  supporting  the  polynomial  complexity  result  for  our  algorithm.  Instead,  top 
scales  linearly  only  for  smaller  problems  (P  =  2),  while  it  has  a  higher-order  polynomial 
behavior  for  larger  problems  (P  =  4)  (the  line  is  still  a  straight  line  but  with  steeper  slope, 
so  it  is  polynomial  on  logarithmic  axes).  This  behavior  depends  on  the  performance  of  the 
chosen  numerical  solver,  and  it  can  improve  benefiting  of  future  advancements  in  the  solver 
implementation.  In  Table  3,  we  compare  the  CP  procedure  with  two  tools,  PRISM  [5]  and 
PARAM  [17],  in  terms  of  runtime,  for  varying  values  of  P  and  R.  Although  neither  tool  solves 
the  same  problem  addressed  in  this  paper,  the  comparison  is  useful  to  assess  the  practicality 
of  the  proposed  approach.  In  particular,  PRISM  only  verifies  PCTL  properties  of  MDPs 
with  no  uncertainties.  PARAM  instead  derives  a  symbolic  expression  of  the  satisfaction 
probabilities  as  a  function  of  the  model  parameters,  to  then  find  the  parameter  values  that 
satisfy  the  property.  Hence,  PRISM  only  considers  a  special  case  of  the  models  considered  in 
this  paper,  while  our  approach  only  returns  the  worst-case  scenario  computed  by  PARAM. 
Results  show  that  the  CP  procedure  runs  faster  than  PRISM  for  some  benchmarks,  but  it 
is  slower  for  larger  models.  This  is  expected  since  the  scalability  of  our  approach  depends 
mainly  on  the  problem  size,  while  the  performance  of  the  iterative  engine  in  PRISM  depends 
on  the  problem  size  and  on  the  number  of  iterations  required  to  achieve  convergence,  which 
is  dependent  on  the  problem  data.  Finally,  our  approach  is  orders  of  magnitude  faster  than 
PARAM,  so  it  should  be  preferred  to  perform  worst-case  analysis  of  system  performances. 

6.2  ZeroConf  Dynamic  Configuration  Protocol  for  IPv4  Link-Local  Addresses 

The  ZeroConf  protocol  [29,30]  is  an  Internet  Protocol  (IP)-based  configuration  protocol  for 
local  (e.g.  domestic)  networks.  In  such  a  local  context,  each  device  should  configure  its  own 
unique  IP  address  when  it  gets  connected  to  the  network,  with  no  user  intervention.  The 
protocol  thus  offers  a  distributed  ’’plug-and-play”  solution  in  which  address  configuration 
is  managed  by  individual  devices  when  they  are  connected  to  the  network.  The  network  is 


Table  3:  Runtime  Comparison 


Tool 

P  =  2,  R  =  2 
N  +  T  =  764 

R  =  7 
2,604 

R  =  128 
47, 132 

P  =  4,  R  =  2 
97, 888 

R  =  32 
1,262,688 

R  =  44 
1,979,488 

P  =  6,R  =  4 
14,211,904 

CP 

0.02s 

0.1s 

2.1s 

8.3s 

1,341s 

2,689 

TO 

PRISM 

0.01s 

0.09s 

196s 

Is 

2, 047s 

TO 

1860s 

PARAM 

22.8s 

657s 

TO 

TO 

TO 

TO 

TO 

composed  of  DVtot  devices.  After  being  connected,  a  new  device  chooses  randomly  an  IP 
address  from  a  pool  of  IPa  —  65024  available  ones,  as  specified  by  the  standard.  The  address 
is  non-utilized  with  probability  po  =  1  —  DVtot / IPa-  It  then  sends  messages  to  the  other 
devices  in  the  network,  asking  whether  the  chosen  IP  address  is  already  in  use.  If  no  reply 
is  received,  the  device  starts  using  the  IP  address,  otherwise  the  process  is  repeated. 

The  protocol  is  both  probabilistic  and  timed:  probability  is  used  in  the  randomized 
selection  of  an  IP  address  and  to  model  the  eventuality  of  message  loss;  timing  defines 
intervals  that  elapse  between  message  retransmissions.  In  [30],  the  protocol  has  been  modeled 
as  an  MDP  using  the  digital  clock  semantic  of  time.  In  this  semantic,  time  is  discretized  in 
a  finite  set  of  epochs  which  are  mapped  to  a  finite  number  of  states  in  an  MDP,  indexed  by 
the  epoch  variable  te.  To  enhance  the  user  experience  and,  in  battery-powered  devices,  to 
save  energy,  it  is  important  to  guarantee  that  a  newly-connected  device  manages  to  select  a 
unique  IP  address  within  a  given  deadline  dl.  For  numerical  reasons,  we  study  the  maximum 
probability  of  not  being  able  to  select  a  valid  address  within  dl.  In  PCTL  syntax: 

P^ax  [ip]  :=  P™ax  (->{ unique  .address }  U  {te  >  dl})  (14) 

We  analyzed  how  network  performances  vary  when  there  is  uncertainty  in  estimating: 
1)  the  probability  of  selecting  an  IP  address,  and;  2)  the  probability  of  message  loss  during 
transmission.  The  former  may  be  biased  in  a  faulty  or  malicious  device.  The  latter  is  esti¬ 
mated  from  empirical  data,  so  it  is  approximated.  Further,  the  IMDP  semantic  of  IDTMCs 
(Section  1),  which  allows  a  nature  to  select  a  different  transition  distribution  at  each  execu¬ 
tion  step,  properly  models  the  time-varying  characteristics  of  the  transmission  channel. 

In  Figure  4,  we  added  uncertainty  only  to  the  probability  of  message  loss  using  the 
likelihood  model,  which  is  suitable  for  empirically-estimated  probabilities.  Using  classical 
results  from  statistics  [11],  we  computed  the  value  of  parameter  (3  from  Set  (2)  corresponding 
to  several  confidence  levels  Cl  in  the  measurements.  In  particular,  0  <  Cl  <  1  and  Cl  = 
1  —  cdfx 2  (2  *  ( Pmax  —/?)),  where  cdfx 2  is  the  cumulative  density  function  of  the  Chi-squared 
distribution  with  d  degrees  of  freedom  (d  =  2  here  because  there  are  two  possible  outcomes, 
message  lost  or  received).  Results  show  that  the  value  of  Pffax  increases  by  up  to  ~10x 
for  decreasing  Cl,  while  classical  model-checking  would  only  report  the  value  for  C'l  =  1, 
which  roughly  over-estimates  network  performance.  The  plot  can  be  used  by  a  designer  to 
choose  dl  to  make  the  protocol  robust  to  varying  channel  conditions,  or  by  a  held  engineer 
to  assess  when  the  collected  measurements  are  enough  to  estimate  network  performances. 

In  Figure  5,  we  compose  different  models  of  uncertainty,  i.e. ,  we  also  add  uncertainty  in 
the  probability  of  selecting  the  new  IP  address  using  the  interval  model.  This  probability  thus 
lies  in  the  interval  [(1  —  u)po,  (1 +  u)po]-  We,  arbitrarily,  fixed  dl  =  25  and  swept  DVtot  in  the 
range  [10—  100],  which  covers  most  domestic  applications,  to  study  how  network  congestion 
affects  the  value  of  Equation  14.  We  studied  four  scenarios:  the  ideal  scenario,  returned  by 
classical  model-checking  techniques;  the  confident,  normal,  conservative  scenarios,  where  we 
added  increasing  uncertainty  to  model  different  knowledge  levels  of  the  network  behavior,  a 


Fig.  4:  Value  of  Equation  14  (top)  and  Fig.  5:  Value  of  Eq.  14  for  increasing 

verification  runtime  (bottom).  number  of  devices  in  the  network. 


situation  that  often  arises  during  the  different  design  phases,  from  conception  to  deployment. 
Results  show  that  P^ax  [tp]  gets  up  to  ~  15  x  higher  than  the  ideal  scenario,  an  information 
that  designers  can  use  to  determine  the  most  sensitive  parameters  of  the  system  and  to  assess 
the  impact  of  their  modeling  assumptions  on  the  estimation  of  network  performances. 

7  Conclusions  and  Future  Work 

We  addressed  the  problem  of  verifying  PCTL  properties  of  Convex-MDPs  (CMDPs),  i.e., 
MDPs  whose  state  transition  probabilities  are  only  known  to  lie  within  convex  uncertainty 
sets.  Using  results  on  strong  duality  for  convex  programs,  we  proved  that  model  checking 
is  decidable  in  PTIME  for  the  fragment  of  PCTL  without  the  Bounded  Until  operator.  For 
the  entire  PCTL  syntax,  the  algorithmic  complexity  becomes  pseudo-polynomial  in  the  size 
of  the  property.  Verification  results  on  two  case  studies  show  that  uncertainty  substantially 
alters  the  computed  probabilities,  thus  revealing  the  importance  of  the  proposed  analysis. 

As  future  work,  we  aim  to  relax  the  rectangular  uncertainty  assumption,  to  limit  the 
power  of  nature  and  obtain  a  less  conservative  analysis.  Also,  we  plan  to  verify  a  complex 
physical  system,  e.g.  an  airplane  power  system,  in  which  modeling  uncertainties  are  present 
both  in  the  underlying  physical  process  and  in  the  failure  probabilities  of  its  components. 
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Appendices 

A  Convex  Optimization  Results 

In  this  appendix,  we  give  details  on  the  results  from  convex  theory  and  duality  that  we  used 
in  the  paper.  Most  of  the  material  is  an  elaboration  of  [11]. 

We  begin  by  giving  the  definition  of  a  convex  set. 

Definition  A.l.  A  set  C  is  convex  if  the  line  segment  between  any  two  points  in  C  lies  in 
C,  i.e.,  if  for  any  x,y  G  C  and  any  a  with  0  <  a  <  1,  we  have:  ax  +  (1  —  a)y  G  C  [24]. 

The  convex  sets  A“,Vs  G  S,\/a  G  -4(s)  model  the  uncertainty  in  the  estimation  of  the  rows 
in  the  transition  matrices  of  A 4c-  In  the  following,  we  will  also  use: 

Definition  A. 2.  A  function  h  :  RN  — >  K.  is  convex  if  its  domain  D  is  a  convex  set,  and  for 
all  x,  y  G  D  and  a  with  0  <  a  <  1,  we  have:  h  (ax  +  (1  —  a)y)  <  ah(x)  +  (1  —  a)h(y)  [24]- 

We  now  introduce  the  convex  uncertainty  models  explicitly  supported  by  our  framework. 
In  order: 

1.  Interval 

2.  Likelihood 

3.  Ellipsoidal 

4.  Entropy  (not  introduced  in  the  paper  for  space  limits) 

For  each  of  them,  we  also  give  details  on  the  primal  and  dual  formulation  of  the  inner 
problem: 

max  xTfsa  (15) 

and  derive  the  time-complexity  of  the  algorithms  used  to  solve  it.  Finally,  for  both  the 
interval  and  ellipsoidal  models  of  uncertainty,  we  provide  the  full  formulation  of  Problem 
(10)  used  to  verify  the  IA  operator  in  polynomial  time.  Similar  results  can  be  obtained  also  for 
the  minimization  problem.  In  the  following,  we  omit  state  and  action  indices  when  possible 
to  improve  readability. 

A.l  Interval  Model 

A  common  description  of  uncertainty  for  transition  matrices  is  by  intervals: 

Ta  =  {fa  g  |  0  <  f“  <  f“  <  f“  <  1,  T'C  =  1}  (16) 

where  f“,fg  G  R*  are  vectors  containing  lower  and  upper  bounds  of  the  elements  of  f.  This 
representation  is  suitable  when  the  components  of  the  transition  matrices  are  individually 
estimated  by  statistical  data. 

We  rewrite  the  inner  problem  in  Equation  (15)  in  primal  form: 

a*  (x)  =  max  xTf 

s.t.  lTf  =  l  (17) 

f  <f  <f 

The  dual  problem  reads: 

cf(x)  =  min  Ai  —  fTA2  +  fTA3 

Al  ,A2 ,A3 

s.t.  A2  A  0,  A3  A  0  (18) 

x  +  A2  —  A3  —  Ail  =  0 

Since  the  primal  problem  is  an  LP,  strong  duality  holds  and  a*  =  d*  [24]. 


Replacing  Problem  (17)  with  Problem  (18),  we  obtain  a  new  LP  formulation  for  Problem 
(10),  used  to  verify  the  U  operator: 

min  xT  1 

X.  A“  s  , A^  s,Ag  s 

s.t.  xs  =  0  Vs  G  Sno 

xs  =  l  Vs  G  Syes  (19) 

-  (f:)TAl,s  +  (TJ'AS,.  Vs  G  S?,  Va  G  A(s) 

x  +  A|  s  -  A|  s  -  A“  S1  =  0  Vs  G  S?,  Va  G  A(s) 

A*  s  >  0,  A|  s  >0  Vs  G  5?,  Va  G  A(s) 

During  the  verification  of  the  Next  operator  instead,  we  want  to  solve  Problem  (18).  To 
derive  the  time  complexity  of  this  operation,  we  rewrite  the  problem  as  [11]: 

d*  =  min  (f  -  f)T  (A1  -  x)+  +  xTf  +  A  (1  -  lrf) 

where  v+  represent  the  positive  part  of  vector  v.  In  this  form,  the  dual  problem  is  uncon¬ 
strained,  and  it  minimizes  a  convex  piecewise  function  with  break-points  at  the  origin  and 
at  Xi,  i  =  1,  •  •  •  N.  A  bisection  algorithm  over  the  discrete  set  b  =  0,  Xi,  i  =  1,  •  •  •  N  will 
thus  find  the  optimal  solution  in  O  ( Nlog{N ))  steps. 

A. 2  Likelihood  Model 

The  likelihood  model  is  appropriate  when  the  transition  probabilities  between  states  are 
determined  experimentally.  The  resulting  empirical  frequencies  of  transition  associated  to 
action  a  €  A  are  collected  in  matrix  Ha.  Uncertainty  in  the  transition  matrices  can  then  be 
described  by  the  likelihood  region  [13]: 

Ta  =  {Fa  G  RWXN  [FVO,  Fa  1  =  1,  £S)S,  hass,  log  (/“ ,)  >  /?“} 

where  fda  <  f3 max  =  £s  s,  h^s,  log (h®s,)  represents  the  uncertainty  level.  Since  the  likelihood 
region  above  does  not  satisfy  Assumption  2.1,  it  must  be  approximated  by  projection  onto 
each  row  of  the  transition  matrix.  We  obtain: 

T:  =  {f“  G  I  f“  >  0,  lTfs“  =  1,  £s,  Ks,  log (/“ ,)  >  pn  (20) 

Even  with  this  approximation,  likelihood  regions  are  less  conservative  uncertainty  represen¬ 
tations  than  intervals,  which  arise  from  further  projection  onto  the  row  components. 

We  rewrite  the  inner  problem  in  Equation  (15)  in  primal  form: 

cr*(x)  =  max  xTfs 

s.t.  lTfs  =  1  (21) 

£s,  hss ,  log (fsa')  >  ps 

fs  >  0 

The  dual  problem  reads  [11]: 

d *(x)  =  min  Ai  -  (1  +  ^S)A2+A2  £s,  hss>  log  (  ) 


s.t.  Ai  >  xmax  =  max  xs 1 
s'es 

A2  >  0 


(22) 


The  primal  problem  is  convex,  and  it  satisfies  Slater’s  condition  [24]  for  non-trivial  un¬ 
certainty  sets,  i.e.  for  (3S  <  /3max  =  E «  s'  ^ss'  log(/is«')i  so  strong  duality  holds  and  a*  =  d* . 
Also,  it  can  be  proven  that  the  dual  Problem  (22)  is  jointly  convex  in  A  and  x. 

Replacing  Problem  (21)  with  Problem  (22),  we  thus  obtain  a  new  formulation  for  Problem 
(10),  used  to  verify  the  U  operator: 

min  x]f  1 

zs,A?jS,A“  „ 


s.t.  xs  =  0 
x„  =  1 


Xu  > 


A  Is  -  (1  +  m,s+^,s  Y,S>  Ks'  log 


A“  s  >  xmax  =  max  xs 

s'GS 

A2,s  >  0 


Vs  G  Sno 

Vs  G  Syes  (23) 
Vs  G  S7  ,V  a  G  A 
Vs  G  S7  ,V  a  G  A 
Vs  G  S7 ,V a  G  A 


Moreover,  when  verifying  the  Next  operator,  the  dual  problem  can  be  reduced  to  one  di¬ 
mension  and  solved  using  a  bisection  algorithm  [11],  with  resulting  time  complexity  O  (Nlog(xmax / e)))  [24] 
with  e  equal  to  the  machine  precision  and  xmax  <  1,  since  x  is  a  vector  of  probabilities. 


A. 3  Ellipsoidal  Model 

Ellipsoidal  models  can  be  seen  as  a  second-order  approximation  of  the  likelihood  model  [11]. 
Intuitively,  in  this  model  the  elements  of  f“  G  J7®  are  restricted  to  lie  on  the  intersection  of 
the  ellipse  E“  =  {f“  |  ||f?“f„a|j2  <  l,i?“  >-  0}  and  the  probability  simplex  AN  =  {f“  G  R"  | 
T'f"  =  1,  f“  >  0}.  We  will  thus  consider  sets: 

Ta  =  {fa  e  RN  |  fa  >  ^  ^fa  =  ^  pa  (fa  _  ha}  ||2  <  1,  0}  (24) 


where  the  matrix  i?“  represents  an  ellipsoidal  approximation  of  the  likelihood  region  r“  = 
{f“  G  R"  |  ]TS,  /i“s,  log(/“s,)  >  /?“}.  In  particular,  can  be  computed  as  follows.  First,  we 
write  the  second  order  approximation  of  the  likelihood  region: 


=  {f“  e 


Es-  Ks'  log(/“S')  {f“  e 


Ey  (°^;s,)2  <  /c2} 


(25) 


where  K?  =  2  (/ 3max  —  /3)  is  a  measure  of  the  uncertainty  in  approximating  the  values  /i“s, . 
The  approximation  in  Equation  25  can  be  written  in  matrix  form: 


Es 


(/>V",)2  <  a:2  ^  IK  (fsQ  -  K)  II2  <  1 


(26) 


R 


a 

s 


(VKsoV-1  0  •••  0 

0  (v^K)-1---  0 

0  •••  : 

0  0 


Matrix  i?“  is  guaranteed  to  be  positive  definite,  i.e.,  i?“  >~  0,  because  it  is  diagonal  and  the 
values  h°ss,  represent  empirical  frequencies  of  transition  from  states  s  to  s' ,  hence  they  are 
non-negative  by  definition. 


We  rewrite  the  inner  problem  in  Equation  (15)  in  primal  form: 


cr*(x)  =  max  xTf“ 

s.t.  lTfs“  =  1  (27) 

IK  (f“  -  K)  ||2  <  1 

ff>o 


The  dual  reads: 


d*  (x)  =  min  Ai  +  A2  +  hTi?.A3 

^1 ,A2,A3 

s.t.  || A3 1| 2  <  A2  (28) 

x  —  Ail  —  i?TA3  =  0 

A2  >  0,  A3  >  0 


where  the  state  and  action  indices  have  been  dropped  to  improve  readability.  The  inner 
problem  is  a  Second  Order  Cone  Problem  (SOCP),  which  satisfies  Slater’s  condition  [24]  for 
any  non-trivial  uncertainty  set,  so  strong  duality  holds  and  a*  =  d* .  The  dual  Problem  (28) 
is  an  SOCP,  so  it  is  trivially  jointly  convex  in  A  and  x. 

Replacing  Problem  (27)  with  Problem  (28),  we  can  thus  obtain  a  new  SOCP  formulation 
for  Problem  (10),  used  to  verify  the  U  operator: 


mm 


s.t.  xs  =  0 

Vs  G  Sno 

xs  =  1 

Vs  G  Sves 

>  A?t.  +  A“  s  +  l«At, 

Vs  G  S'7,  Va 

11^3, s  2  <  ^2 ,s 

Vs  G  S7,  Va 

x  -  A?iSl  -  i?“TA|  s  =  0 

Vs  G  S7,  Va 

A£,s  >  0,  A|  s  >  0 

Vs  G  S7,  Va 

(29) 


Introducing  uncertainty  comes  at  the  cost  of  solving  an  SOCP  with  ( N  +  2 )MQ  more 
variables  and  ( N  +  1)MQ  more  constraints  than  the  original  LP,  where  Q  =  |S7|  =  O(N). 

During  the  verification  of  the  Next  operator  instead,  we  want  to  solve  either  Problem  (27) 
or  Problem  (28).  Since  they  are  both  SOCP,  they  can  be  solved  using  interior-point  methods 
with  worst-case  ( practical )  time  complexity  O  ( N1'5log(xmax  \  e))  ( O  ( Nlog(xmax  \  e)))  [24] 
and  xmax  <  1,  since  x  is  a  vector  of  probabilities. 

Finally,  we  report  a  second  dual  formulation  of  Equation  (15),  which  we  experimentally 
found  to  have  better  runtime  performance  when  using  the  VI  routine  to  verify  properties 
containing  the  IA  operator.  Intuitively,  the  faster  performance  is  achieved  because  this  for¬ 
mulation  allows  to  write  the  analytical  expression  of  the  dual  solution  of  the  inner  problem 
as  a  function  of  the  state  probabilities  x.  The  analytical  expression  can  then  be  rapidly  eval¬ 
uated  during  the  VI  iterations  using  the  values  of  x  estimated  in  the  previous  iteration.  We 
note  that  the  same  approach  cannot  be  used  for  the  CP  routine,  since  the  derived  analytical 
expression  is  not  easily  representable  in  conic  form  when  the  optimization  Problem  (29) 
operates  both  on  the  decision  variables  x  and  A  at  the  same  time.  In  the  following,  we  will 
drop  the  state  and  action  indices  to  improve  readability. 


We  rewrite  the  inner  problem  in  Equation  (27)  in  an  equivalent  primal  form: 


a*  (x)  =  max  xTf 
f 

s.t.  lrf  =  1 


ifs'-h,'? 


EUs'  —  »-s' 


<  1C 2 


f  >  0 


The  Lagrangian  operator  associated  to  Problem  (30)  reads: 


£(f,/b£,  v)  =  xTf  +  A‘(l  -  lTf)  +^Tf  +  ^  ^ S>  h  hs^ 


(30) 


(31) 


The  primal  optimal  value  can  be  computed  as: 

a*  (x)  =  max  min  C{ f ,  fi,  £,  v)  (32) 

f  n,i,v 

By  the  minimax  theorem  [24],  we  obtain  an  upper  bound  on  the  value  of  cr*(x)  by 
inverting  the  “max”  and  “min”  operators: 


d* (x)  =  min  max  C( f ,  /z,  £,  v) 


(33) 


Problem  (30)  satisfies  Slater’s  condition  [24]  for  any  non-trivial  uncertainty  set,  so  strong 
duality  holds  and  a*  =  d*.  We  can  thus  solve  Problem  (33)  instead  of  Problem  (30)  while 
preserving  the  soundness  and  completeness  of  the  verification  procedure.  As  a  first  step,  we 
compute  the  dual  function  g(n,  £,  v)  by  solving  the  inner  problem  in  Equation  (33),  i.e. ,  we 
aim  at  computing: 

£.,”)  =  max£(f)  (34) 

We  can  solve  Problem  (34)  by  setting  the  gradient  of  the  Lagrangian  to  zero,  and  solving 
for  the  optimal  primal  solution  f*: 
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Substituting  f*  back  into  Problem  34,  we  obtain  the  dual  function: 

0,z/>  0)  =  ^  +  tW2  +  ^](/v(av  -M +  £«'))  +  ^^2(hs'(xs'  ~  g  +  L')2)  (36) 


We  can  now  compute  d*  solving  the  dual  problem: 


d*  =  min  g{n,£,u) 


(37) 


r 


Before  further  proceeding,  we  note  that,  for  nronotonicity  reasons,  we  can  trivially  set 
=  0.  We  thus  aim  to  solve  the  following  optimization  problem: 

d*  =  min^  fj,  +  uK2  +  ^ ](hs>(xs >  -  /x))  +  —  y ^(h3>(xa>  -  g)2)  (38) 

—  s'  s' 


We  first  set  the  partial  derivatives  of  the  dual  function  g  to  zero  and  compute  the  optimal 
dual  solution  (g*,i/*).  Formally: 


U  =  l-g(*-Kl)T-Vh  =  0 


H*  =  Es'  hs'Xs' 

V*  =  2^a/Es'  M^V  -  g)2 


(39) 


The  optimal  value  can  then  be  computed  as  d*  =  g{g* ,  v*). 

Finally,  we  report  the  primal  and  dual  formulations  to  solve  the  inner  problem  when 
verifying  properties  of  the  form  P>p  [tp\.  In  the  primal  problem  of  Equation  (15),  we  change 
the  optimization  operator  from  “max”  to  “min” : 


<t*  (x)  =  min  xTf 
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s.t.  lTfsa  =  1 
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f  >  0 


(40) 


Following  the  same  steps  presented  above,  we  obtain  the  corresponding  dual  problem: 


d*  =  max  g  —  vK' 


+  X](Mav  -  m))  -  7 V-  EtMav  -  v)2) 


which  admits  the  optimal  solution: 

fi  =  1+!b(x-M)T-lrh  =  0 

I  §9.  _  ,  h(x-^l)2T  _  n 

l  di>  ~  ^  +  4v2  ~  U 


M*  =  Es'  hs'Xs' 

V*  =  aWEs'  hs>(xa>  -  g)2 


(41) 


(42) 


A. 4  Entropy  Model 

The  entropy  model  of  uncertainty  can  be  viewed  as  a  variation  of  the  likelihood  model. 
In  the  likelihood  setting  we  bound  the  divergence  from  an  empirically  extracted  distribu¬ 
tion,  whereas  in  the  entropy  setting  we  bound  the  divergence  from  a  reference  analytical 
distribution  q  [11].  We  will  thus  consider  sets: 

T:  =  {f“  G  r  I  f“  >  0,  lrf“  =  1,  ES'  fss'  log  (§;)  <  /?“}  (43) 

We  rewrite  the  inner  problem  in  Equation  (15)  in  primal  form: 

*  T  X* 

(T  =  max  x  r  s 

s.t.  lTfs  =  1 

E«'  fss'  log  (  ^)  <  Pa 

is  >  0 


(44) 


The  dual  problem  reads: 


d*  =  mm  A  log  (X)g,  qss>  exp  (^))  +  /3SA  (45) 

s.t.  A  >  0 

The  primal  problem  is  convex,  and  it  satisfies  Slater’s  condition  [24]  for  non-trivial  un¬ 
certainty  sets,  i.e.  for  /3S  >  0,  so  strong  duality  holds  and  cr*  =  d*. 

Replacing  Problem  (44)  with  Problem  (45),  we  can  thus  obtain  a  new  formulation  for 
Problem  (10),  used  to  verify  the  U  operator: 
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We  prove  its  joint  convexity  in  x  and  A®  as  follows.  The  cost  function  and  Constraints  (46a), 
(46b)  and  (46d)  are  trivially  convex.  Constraint  (46c)  is  generated  by  a  primal-dual  transfor¬ 
mation,  so,  according  to  convex  theory,  it  is  convex  in  the  dual  variables  A®  by  construction. 
Moreover,  we  prove  that  it  is  also  jointly  convex  in  x  by  induction  on  the  number  TVS'  of 
next  states  s'  G  S  for  state  s.  As  a  base  case,  TVS  =  1,  and  we  can  rewrite  the  constraint  as: 

xs  >  A®  log  (q8a>)  +  xs‘  +  /3SA“  (47) 

which  is  trivially  jointly  convex.  We  now  assume  that  the  constraint  is  jointly  convex  for 
N S  =  n,  and  prove  that  it  is  jointly  convex  also  for  NS  =  n  +  1.  This  result  immediately 
follows  from  observing  that  increasing  NS  simply  introduces  one  more  term  in  the  summa¬ 
tion,  so  if  the  constraint  is  jointly  convex  for  TVS'  =  n  it  must  remain  jointly  convex  also 
for  TVS  =  71  +  1,  since  an  affine  addition  preserves  convexity  according  to  convex  theory.  We 
conclude  that  Problem  (46)  is  convex. 

Moreover,  when  verifying  the  Next  operator,  the  dual  problem  is  unidimensional  and 
it  can  thus  be  solved  using  a  bisection  algorithm  [11],  with  resulting  time  complexity 
O  (Nlog(xmax\  e)))  [24]  with  e  equal  to  the  machine  precision  and  xmax  <  1,  since  x  is 
a  vector  of  probabilities. 


B  Proof  of  Contraction  Lemma 


In  this  appendix,  we  present  a  verification  procedure  for  the  Until  operator  based  on  Value 
Iteration  (VI),  which  should  be  considered  as  an  alternative  to  the  CP  procedure.  We  report 
also  the  VI  procedure  because  it  has  shorter  runtime  than  the  CP  procedure  for  some  of  the 
analyzed  case  studies,  depending  on  the  structure  and  data  of  the  problem.  Both  procedures 
should  be  run  in  the  early  stages  of  system  verification  and  the  one  which  performs  better 
should  be  used  in  the  rest  of  the  project  development. 

We  start  by  defining: 

Definition  B.l.  Contraction.  Let  ( B,d )  be  a  metric  space  and  g  :  B  — ►  B.  Function  g  is  a 
contraction  if  there  is  a  real  number  9,  0  <  8  <  1,  such  that: 

d  (g(u),  g(v))  <  9d(u,v)  \/u,v  G  B  (48) 


In  the  following,  we  will  use: 

Proposition  B.l.  Contraction  mapping.  Let  ( B ,  d)  be  a  complete  metric  space  and  g  :  B  — ► 
B  a  contraction.  Then  there  exists  a  unique  point  x*  G  B  such  that: 


g(x*)  =  x* 


Additionally,  if  x  £  B,  then: 

lim  gk( x)  =  x* 

/c— >-+oo 

We  use  the  mapping  g  =  G  defined  as: 


0  Vs  G  Sno 


G  = 


1  Vs  G  Syes 


0  Vs  G  S'  A  i  =  0 

max  max  (x®-1)Tff  Vs  G  S'  A  i  >  0 
aeA(s)ff-e^ 


(49) 


where  Syes  d=f  Sat  (P>i  [faUfa]),  Sno  d=f  Sat  (P^iU^])  and  S?  =  S  \  (, Sno  U  S'"68).  We 
note  that  MQ  convex  problems  need  to  be  solved  to  compute  mapping  G,  with  Q  =\  S?  |.  For 
the  uncertainty  models  considered  in  the  paper,  each  problem  can  be  solved  with  complexity 
O  (Vlog(l/e))  [11],  and  e  equal  to  the  machine  precision.  To  simplify  notation  in  the  proof, 
we  use  the  weighted  maximum  norm  ||  .  ||w  of  a  vector  v  G  defined  as: 


INI 


m  ax  - 

i=l. ..JV  Wi 


where  Wi  is  the  scalar  weight  associated  to  each  element  of  v. 
We  can  now  state: 


(50) 


Lemma  B.l.  Mapping  G  is  a  contraction  over  the  metric  space  (R",  ||  .  ||w). 

Proof.  The  proof  follows  closely  the  ones  in  [B.l]  (Vol.  II,  Section  2.4)  and  [20].  Those  proofs 
refer  to  a  control  setting,  where  the  optimal  action  (control)  can  be  selected.  Hence,  the 
contraction  needs  to  hold  for  only  one  of  the  available  actions,  i.e.  the  optimal  one  (existential 
quantification).  Conversely,  in  the  verification  setting,  the  contraction  needs  to  hold  across  all 
available  actions,  because  we  consider  the  worst  case  resolution  of  nondeterminism  (universal 


quantification).  Further,  as  in  [20],  we  quantify  across  all  nature  behaviors:  this  is  possible 
due  to  Assumption  2.2.  For  the  sake  of  brevity,  in  the  following  we  will  only  consider  the 
calculation  of  Pr™m,  but  the  same  reasoning  applies  also  for  the  maximization  problem. 

We  start  from  partitioning  the  state  space  S  =  Syes  U  Sno  U  S?  as  explained  in  Section 
5.2.  Since  at  all  iterations  the  probabilities  Pr™*"  will  remain  constant  by  construction  in 
all  states  s  G  Syes  U  Sno,  we  do  not  need  to  consider  these  states  explicitly.  In  particular, 
we  perform  the  following  transformations  of  the  CMDP  underlying  graph:  we  collapse  the 
set  Syes  into  one  terminal  state  t,  and  eliminate  all  states  s  €  Sno  from  the  graph.  These 
transformations  are  fundamental  together  with  Assumption  2.2  to  guarantee  that  all  possible 
adversaries  Adv  are  proper  in  the  transformed  graph,  i.e.  they  almost  surely  reach  the 
terminal  state  t  for  all  transition  matrices  in  T  [31].  We  will  now  work  with  the  new  state 
space  =  S’-  U  {t},  and,  for  simplicity,  we  redefine  N  =  \  S t  |.  We  further  partition  ,  as 
follows.  Let  Si  =  {t}  and  for  q  =  2, 3,  •  •  •  compute: 

S„  =  {s  €  |  s  4.  Si  U  •  •  •  U  S„_i,  min  max  min  /“  <  >  0) 

91  1  ^  9  1  aeA(s)  s'es1u-usq-1f-e^fss  1 


Let  r  be  the  largest  integer  such  that  Sr  is  nonempty.  Since  all  adversaries  are  proper,  we 
are  guaranteed  that  U q—iSq  =  SA  We  now  need  to  choose  weights  ws,\/s  €  S*  such  that  G 
is  a  contraction  with  respect  to  ||.||w.  First,  we  take  the  sth  component  ws  to  be  the  same 
for  states  s  in  the  same  set  Sq.  Then  we  set  ws  =  yq  if  i  €  Sq,  where  y-\ ,  •  •  •  ,yr  are  scalars 
satisfying  1  =  y\  <  y-2  <  ■  ■  ■  <  yr-  Further,  let: 
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mm  mm  mm  mm 
q= 2,---  ,r  a£A  s€Sq  ff'GJ7" 
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s'eSiU-us,-! 


By  construction  0  <  £  <  1. 

The  rest  of  the  proof  goes  as  follows:  first,  we  will  show  that  if  we  can  find  y-2 ,  •  ■  ■  ,  yr 
such  that  for  q  =  2,  •  •  ■  ,  r: 

<Q 

Vq  Vq 

for  some  6  <  1,  then  G  is  a  contraction.  Second,  we  will  prove  that  such  values  always  exist. 
We  begin  by  defining: 
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G?(x) 


min  min  xTf“ 

aGyt(s)  f“e ^ 

min  xTf“ 

f»e^“ 


i.e.  the  sth  element  of  the  output  of  mapping  G  applied  to  vector  x  G  RN,  and  the  same 
element  when  mapping  G  is  evaluated  only  at  the  fixed  action  a  G  A(s).  Then,  for  all  vectors 
v,  u  G  RN,  we  determine  A(s)  such  that: 

a  =  argmin  G{u ) 

Ms) 


We  can  thus  write  for  all  sG5*: 

G»(v)-Gs(u)  =  G8(v)-G“(u) 
<G“(v)-G“(  u) 

=  ES'  (V>.'  -  u:s,us,) 

<  E S'  Mss‘  -  US') 


where: 


V“  =  argmin  vTf“ 

f“G  J7? 

U“  =  argmin  uTf“ 

Mss’  =  argmax{V“s,  (tv  -  us>) ,  U“s,  (vs,  -  tv)} 


Let  q(s)  be  such  that  state  s  belongs  to  the  set  Sq(sy  Then,  for  any  constant  c: 

II v  -  u||w  ^vs-us<  cyq(s)  Ms  G  S 1 


We  can  thus  write  Ms  G  Sq  and  q  =  1,  •  •  •  ,  r: 
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We  have  thus  proved  that  G*(v)  G*(u)  <  c(^  for  an  arbitrary  state  s  G  Si  Taking  the 
maximum  over  St,  we  get: 


|| G'(tt)  —  G(tt)||w  <  c$,  Vu,v  G  R"s.t.||v  —  u||  <  c 


so,  G  is  a  contraction  over  the  metric  space  (KN,  ||  .  ||w),  and: 


6=  max  ^(l-0  +  —  (51) 

i<9<ryq 

is  the  corresponding  contraction  factor.  Finally,  we  constructively  prove  by  induction  that 
it  is  always  possible  to  find  scalars  y i,  •  •  •  ,  yr  such  that  the  above  assumptions  hold.  As  the 
base  case,  we  set  yo  =  0,yi  =  1.  At  the  induction  step,  we  suppose  that  y-2 ,  •  •  ■  ,yq  have 
already  been  determined.  If  £  =  1,  we  set  yq+i  =  yq  +  1.  If  £  <1,  we  set  yq+\  =  |  {yq  +  mq) 
where: 

mq  =  i^-<9  {yi  +  _  y*-1)} 

With  these  choices,  we  are  guaranteed  that: 


mq+i  =  min 


|  mq,yi  + 


j 

1-C 


(y* 


y*-i) 


so  by  induction,  we  have  that  yq  <  yq+i  <  mg+i,  and  we  can  construct  the  required  sequence. 

□ 


We  now  state  the  main  results  of  this  appendix: 

Lemma  B.2.  The  VI  procedure  to  verify  the  Unbounded  Until  operator  is  sound  and  com¬ 
plete,  i.e., 

Vmax[(j)1U(j)2}=  lim  Gfc(x) 

k — ^-|-oo 


(52) 


In  practice,  we  need  a  criterion  to  stop  the  iteration,  so  an  error  is  introduced  in  Equa¬ 
tion  (52)  and  Pmax  is  approximated.  In  the  following,  we  show  how  to  compute  an  exact 
lower  bound  Kmm  on  the  number  of  iterations  required  to  obtain  the  desired  accuracy  ed- 
The  existence  of  such  a  lower  bound  also  proves  the  soundness  and  completeness  of  the  VI 
procedure.  Further,  we  propose  a  heuristic  stopping  criterion  based  on  a  relative  tolerance 
check,  which  is  often  used  to  reduce  runtime  [5].  The  user  can  choose  between  the  two 
approaches,  depending  on  the  application. 


Exact  Lower  Bound.  From  Lemma  B.l,  at  the  end  of  the  ith  iteration  the  residual  error 
in  estimation  is  bounded  by: 


Pi  =  HP”10®  -x‘||w  < 

where  po  is  the  initial  error  in  estimation  which  can  be  trivially  bounded  by  po  <  1.  The 
error  in  the  estimation  of  the  satisfaction  probabilities  is  bounded  by  ed  if: 

Pi  <  =H  PTX  -  |<  ed,  Vs  G  S7 

W max 

where  wmax  is  the  maximum  of  the  weights  of  Norm  (50) .  We  can  thus  obtain  a  lower  bound 
Krnln  for  the  number  of  iterations  required  to  achieve  the  desired  accuracy  ed- 

°K  /  .Ts^Tsmin  log[ed(l  -  9)]  -  log{wmax) 

i - a  -  - >  K  >  K  =  -  -  (5 3) 

1  0  Wmax  logi^O) 

Heuristic  Based  on  Relative  Tolerance.  Although  provably  exact,  the  lower  bound 
AT7™”  jn  the  number  of  iterations  derived  in  the  previous  section  might  be  too  conservative 
and  result  in  an  unnecessary  increase  in  runtime.  We  thus  also  present  a  heuristic  stopping 
criterion  based  on  relative  tolerance.  In  particular,  we  stop  when: 

Sr  >  max  (\x\  -  xl~1\/x1^)  (54) 

sGo  • 

the  maximum  relative  difference  in  the  computed  value  of  P'™ax[(t>iU(t>2 ],  Vs  G  S'  between 
two  consecutive  iterations  is  below  a  user  defined  tolerance  5rei .  We  note  that  Such  a  criterion 
does  not  guarantee  that  the  error  is  bounded  by  Sr.  The  required  Sr  to  achieve  accuracy 
erj,  i.e. ,  S*,  depends  on  the  CMDP  model  and  needs  to  be  determined  by  trial-and-error, 
a  common  practice  in  iterative  procedures  (e.g.  the  ODE  solver  in  a  circuit  simulator).  To 
determine  S*,  we  compute  several  approximations  of  Pmax[i/;]  while  decreasing  Sr  by  steps  of 
10 x.  We  heuristically  stop  when  no  probability  Psma;c,Vs  G  S  .  changes  more  than  after 
checking  Criterion  (54)  for  5*  and  (5*/100.  Finally,  errors  in  solving  the  inner  problems,  as 
introduced  in  Section  5.3,  are  propagated  across  iterations.  We  call  einn  the  inner  problem 
accuracy.  If  the  VI  procedure  exits  after  /  iterations  and  ed  <  I  x  einn,  the  procedure  needs 
to  be  run  again  after  decreasing  elnn  to,  approximately,  €inn  <  ed/I- 

We  use  the  VI  routine  with  ed  =  10-3  to  verify  again  <j>  =  P>o.3[  DU  oj]  in  the  example 
in  Figure  1.  After  3  iterations,  we  get  Pmm[  D  U  uj  }  =  [0.2,  0, 1, 0.32]  and  Sat{(j> )  =  {s2,  S3}. 


[B.1]  D.  Bertsekas,  “Dynamic  Programming  and  Optimal  Control”,  Athena  Scientific,  2011 


C  Toy  LP 


This  appendix  reports  the  full  LP  formulation  that  was  used  to  verify  property  </>  =  P>o.3[  Ww] 
on  the  example  in  Figure  1.  Problem  (10)  written  with  the  data  of  the  model  has  19  variables 
and  11  constraints.  All  variables  are  (implicitly)  bounded  to  be  positive  apart  from  the  ones 
labeled  as  free  at  the  bottom  of  the  formulation. 
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Minimuum  Probability  of  Eating 


D  Dining  Philosophers 

We  analyze  the  classical  Dining  Philosopher  Problem  [D.l].  Briefly,  n  philosophers  are  sitting 
at  a  table  with  n  available  forks.  Each  philosopher  can  either  think  or  eat:  when  he  becomes 
hungry ,  he  needs  to  pick  both  the  fork  on  his  right  and  on  his  left  before  eating.  Since 
there  are  not  enough  forks  to  allow  all  philosophers  to  eat  together,  they  need  to  follow 
steps  according  to  a  stochastic  protocol  to  eat  in  turns.  We  consider  this  case  study  relevant 
because  it  can  be  used  to  model  real  shared-resources  stochastic  protocols  [D.  1],  and  because 
the  size  of  the  model  n  can  be  easily  scaled  to  benchmark  the  time  complexity  of  our  routines. 

We  model  the  uncertainty  of  the  philosophers  in  deciding  which  fork  to  pick  first:  while 
the  nominal  protocol  assigns  0.5  —  0.5  probability  to  the  left  and  right  fork,  we  assume 
that  these  values  are  only  known  with  ±i0%  confidence.  The  parameters  for  each  model  of 
uncertainty  corresponding  to  this  level  of  confidence  can  be  set  using  the  approach  suggested 
in  [11].  For  example,  for  the  Interval  model,  the  probabilities  he  in  the  interval  [45%  —  55%]. 
Within  this  setting,  we  aim  to  determine  which  is  the  quantitative  minimum  probability  for 
any  philosopher  to  eat  within  k  steps  of  the  protocol  after  he  becomes  hungry.  In  PCTL 
syntax: 

pmin  ^  ,=  p min  [p  <*  {Eating}]  (56) 

with  initial  states  So  =  Sat(hungry).  Figure  2  shows  the  evolution  of  the  probability  of 
Equation  (56)  as  a  function  of  the  number  of  protocol  steps  k.  As  expected,  the  probability 
of  eating  steadily  increases  as  the  number  of  steps  increases.  However,  the  plot  also  shows 
that  adding  uncertainty  decreases  this  probability  with  respect  to  the  nominal  scenario  (if 
no  uncertainty  is  added,  our  results  match  those  in  [D.l].  The  inset  of  Figure  2  shows  the 
relative  deviation  in  probability  with  respect  to  the  nominal  case:  a  ±10%  uncertainty  can 
cause  a  deviation  up  to  35%  in  the  computed  probabilities,  and  the  deviation  is  always 
higher  than  10%  for  k  <  60.  Further,  the  deviation  is  larger  for  the  Interval  and  Ellipsoidal 
models,  since  they  are  the  most  conservative  among  the  considered  ones,  as  explained  in 
Section  2.2  and  Appendix  A. 

Lastly,  we  evaluate  the  runtime  performance  of  our  routines  while  varying  the  size  n 
of  the  problem.  We  report  three  data  points  for  each  uncertainty  model,  corresponding  to 
n  =  3, 4, 5.  Figure  3  shows  that  runtime  increases  approximately  linearly  with  the  number  of 
states  N.  The  discrepancy  with  the  expected  quadratic  behavior  can  be  explained  considering 
that  in  this  case  study  (and  in  most  practical  ones)  not  all  actions  a  £  A  are  available  at  each 
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Fig.  2:  Evolution  of  Equation  (56)  for  in¬ 
creasing  k  and  n  =  3  philosophers. 


k  =  150  steps. 


state  s  £  S  and  the  transition  matrix  Fa  £  Ta  is  sparse.  The  interval  and  ellipsoidal  models 
run  faster  because  the  inner  convex  optimization  problems  can  be  solved  using  simpler 
atomic  operations  (sum  and  multiplication)  than  the  likelihood  model  (logarithm).  Further, 
the  routine  for  the  interval  model  runs  only  1.2  x  slower  than  the  Certain  Slow  routine, 
and  the  penalty  rises  to  20  x  with  respect  to  the  Certain  Fast  routine:  this  result  can  be 
interpreted  as  the  cost  of  not  being  able  to  use  optimized  library  procedures  for  matrix- vector 
multiplication  when  adding  uncertainty  to  the  model.  These  results  support  our  claim  of 
good  scalability  of  the  proposed  approach  with  respect  to  the  model  size. 
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